- Email: [email protected]
Set-Membership Estimation of Hybrid Systems via SAT Modulo ODE*
16th IFAC Symposium on System Identification The International Federation of Automatic Control Brussels, Belgium. July 11-13, 2012
Set-Membership Estimation of Hybrid Systems via SAT Modulo ODE ? Andreas Eggers ∗ Nacim Ramdani ∗∗ Nedialko S. Nedialkov ∗∗∗ Martin Fr¨ anzle ∗ ∗
Carl von Ossietzky Universit¨ at, Oldenburg, Germany, (e-mail: {eggers—fraenzle}@informatik.uni-oldenburg.de). ∗∗ Universit´e d’Orl´eans, PRISME, 18020 Bourges, France, (e-mail: [email protected]) ∗∗∗ McMaster University, Hamilton, Ontario, Canada, (e-mail: [email protected]) Abstract: Set membership estimation (SME) of nonlinear hybrid systems is still a challenging issue. Although SME of nonlinear continuous systems has made significant progress recently, the direct extension of these methods to the hybrid case is not easy. Meanwhile, satisfiability (SAT) checkers for Boolean combinations of arithmetic constraints over real- and integer-valued variables have made significant progress, as they can effectively deal with algebraic constraints between variables and non-linear ODEs, what is denoted as SAT Modulo ODE. Finally, the corresponding solvers solve in a natural way the hybrid differential and algebraic constraints satisfaction problems that underlie SME of hybrid systems. This paper presents the application of such a SAT Modulo ODE solver to SME of hybrid dynamical systems. Keywords: Bounded error, constraint satisfaction problems, nonlinear systems, uncertain dynamic systems, state estimation, robust estimation, satisfiability. 1. INTRODUCTION State estimation is an important problem when addressing control or diagnosis issues of dynamical systems. Most controlled systems exhibit both smooth continuous dynamics and abrupt switches, hence can be efficiently modeled using hybrid automata, which combine discrete and continuous variables. In practice, the whole hybrid state remains only partially observable, as some switches and some continuous components cannot be directly observed. Hybrid state estimation aims at reconstructing both the discrete mode, hence the switching sequence, and the corresponding continuous state variables, based on a set of possibly discrete measurements, the knowledge of the hybrid model, and assumptions about the uncertainties and perturbations acting on the system. In many cases, it is more natural to assume that all uncertain quantities—measurement noise, model uncertainty and modelling errors—belong to a known set. This assumption is rather natural and requires much less data than statistical assumptions. In such an unknown-but-bounded-error (UBBE) framework, the estimation problem no longer has a unique solution, but there exists a set of state vectors that are consistent with measured data, the model structure and the prior error bounds. Set-membership estimation (SME) techniques allow the characterization of this solution set. SME tech? This work has been supported by the German Research Council DFG within SFB/TR 14 “Automatic Verification and Analysis of Complex Systems” (www.avacs.org), by the French National Research Agency under contract ANR 2011 INS 006 04 “MAGIC-SPS” (projects.laas.fr/ANR-MAGIC-SPS), and by the Natural Sciences and Engineering Research Council of Canada.
978-3-902823-06-9/12/$20.00 © 2012 IFAC
niques have reached a mature stage when dealing with linear or nonlinear discrete-time models; see a.o. (Milanese and Novara, 2011; Kieffer and Walter, 2011; Jaulin, 2011; Le Bars et al., 2012), and the references therein. Bemporad et al. (2005) addressed SME with piecewise affine systems. More recently, SME with continuous-time nonlinear systems has made significant progress; a.o. Moisan et al. (2009); Meslem et al. (2010); Meslem and Ramdani (2011) introduced new enclosing methods for nonlinear ordinary differential equations (ODE) for SME, and Goldsztejn et al. (2010) showed how to embed ODE into the numerical constraint satisfaction problem framework. SME with non linear hybrid systems is a challenging issue. Yet, nonlinear hybrid systems identification using point techniques is prohibitively time consuming, as it usually requires solving combinatorial nonlinear programming or optimization problems. The issue becomes even more difficult in set-membership framework, as SME solving techniques usually need to compute solution space coverings by means of multi-dimensional boxes. To the best of our knowledge, the only work addressing SME with nonlinear hybrid systems is by (Benazera and Trave-Massuyes, 2009), who addressed hybrid systems with discrete-time only nonlinear continuous dynamics. In this paper, we address SME of nonlinear hybrid systems, whose continuous dynamics are described using nonlinear ODEs. We show that SME of hybrid systems boils down to solving differential and algebraic constraints satisfaction problems, which in turn can be written as SAT Modulo ODE formulas, hence solved using appropriate satisfiability checkers— notions to be defined in following sections.
440
10.3182/20120711-3-BE-2027.00292
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
Outline. Sect. 2 recalls hybrid system modeling and SME. Sect. 3 introduces the main contribution of this work, namely the use of satisfiability checkers for solving SME with hybrid systems. Sect. 4 gathers experimental evaluation on two examples, including a hybrid system. 2. SET-MEMBERSHIP ESTIMATION OF HYBRID SYSTEMS We consider a hybrid automaton (Alur et al., 1995) given by H = (Q, D, P, Σ, A, Inv, F). (1) Here Q is a set of locations. Given a location q ∈ Q, the continuous dynamics, and hence flow transitions, are described by non-autonomous differential equations fq ∈ F of the form flow(q) : x(t) ˙ = fq (x, p, t), (2) where fq : D × P × R+ 7→ D is nonlinear and assumed sufficiently smooth over D ⊆ Rn , with dimension n that may depend on q, and p ∈ P, where P is an uncertainty domain for the parameter vector p. Inv is an invariant, which assigns a domain to the continuous state space of each location. It is defined by the following system of inequalities: Inv(q) : νq (x(t), p, t) < 0, (3) where νq : D × P × R+ 7→ Rm is a vector-valued nonlinear function, the negativity constraint applies componentwise, and the number m of inequalities may also depend on q. The set A is the set of discrete transitions {e = (q → q 0 )} each given by the 5-uple (q, guard, σ, ρ, q 0 ), where q and q 0 represent upstream and downstream locations respectively; guard is a condition given as the system of equations guard(e) : γe (x(t), p, t) = 0, (4) σ is an event, and ρ is a reset function. A transition q → q 0 may occur when the continuous flow resides within the guard set, i.e. when the continuous state satisfies condition (4). Let us also consider the following measurement equation y(t) = gq (x(t), p, t), (5) where function gq : D × P × R+ 7→ Rm is nonlinear. We assume that one has access to measurements of the output vector y(t). The latter is usually measured at some sampling time instants tj taken over a time-grid t0 < t1 < t2 < . . . < tnT . The sampling intervals need not be constant, as it may be the case when measuring for instance reactant concentrations in biochemical reactors. The output vector may also be measured using discrete sensors, i.e. sensors that monitor only a single threshold value such as critical fluid level in chemical processes or proximity sensors in robotics applications. Discrete sensors are suitable for monitoring guard conditions that trigger discrete transitions in hybrid systems as emphasized by Koutsoukos (2003). Measurement noise is taken as a discrete time signal assumed additive and bounded with known bounds. Denote by y i the feasible domain for output vector yi at time ti . Recall that the parameter vector p is uncertain and lay within a bounded set with known bounds P. We assume that the initial hybrid state vector is only partially known, i.e. the actual discrete mode may be unknown or the continuous component taken in a bounded domain x(t0 ) ∈ x0 .
The goal of the bounded error estimation is to compute conservative outer enclosures for feasible sets for both the discrete modes and associated continuous variables that are consistent with the feasible domains for measurements, the chosen hybrid model, and the assumptions about the uncertainties and perturbations acting on the system. That is, given two measurements y i and y i+1 gathered at the two time instants ti and ti+1 , the estimation problem boils down to simultaneously: (i) Detect and identify all possible discrete transitions (t? , e) that may occur at time t? ∈ [ti , ti+1 ]. Because there are uncertainties in either initial state vector, continuous dynamics or measurements, there might be a continuum of time instants, where events may occur. There may also exist several such time intervals in [ti , ti+1 ]. (ii) Reconstruct the sets of hybrid states (q(ti ), x(ti )) and (q(ti+1 ), x(ti+1 )) that are consistent with the switching sequence reconstructed in (i), the modelling (1)-(5), and the inclusions y(ti ) ∈ y i and y(ti+1 ) ∈ y i+1 . These two steps must be done for all i, 0 ≤ i ≤ nT − 1. Finally, this set membership estimation problem for hybrid systems can naturally be stated as finding feasible solution sets for Boolean combinations of both arithmetic and differential constraints over real- and integer-valued variables. This is the purpose of the iSAT-ODE solver we introduce in the next section. 3. THE ISAT-ODE SOLVER In (Eggers et al., 2011), we presented an extension of the iSAT solver (Fr¨anzle et al., 2007) with a combination of enclosure methods for ODEs. In this section we concentrate on the application of iSAT-ODE on the set estimation problem and refer the reader to the above references for algorithmic details. 3.1 Encoding and Analysis of Hybrid Systems Our approach to constraint-based analysis of hybrid systems is based on a predicative encoding of the hybrid system’s trajectories. This encoding consists of a predicate describing the valid initial states, a k-fold unwinding of a formula connecting pairs of states by continuous and discrete transitions and—when considering classical bounded model checking problems—a target predicate that describes the final state of the system, whose reachability is to be analyzed. SAT Modulo ODE. Let Φ be a SAT Modulo ODE formula. This is a quantifier-free Boolean combination of arithmetic constraints over bounded real-, integer-, and Booleanvalued variables, simple bounds, and ODE constraints over real variables with the following properties:
441
• arithmetic constraints over variables x, y, and z are of the form x ∼ ◦(y, z) or x ∼ ◦(y), where ∼ is a relational operator from {<, ≤, =, ≥, >}, and ◦ is a total unary or binary operator from {+, −, ·, sin, cos, powN , exp, min, max}; • simple bounds are of the form x ∼ c with ∼ as above a relational operator, x a variable, and c a constant; and
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
• ODE constraints are sufficiently smooth, time invariant, and given by x˙ i = dxi /dt = f (x1 , . . . , xn ) with all occurring variables xi themselves being defined by ODE constraints and f being a function com√ posed of {+, −, ·, /, powN , exp, ln, , sin, cos}. These ODE constraints must occur only under an even number of negations in the formula, allowing e.g. m1 → ((x˙ = sin(y)) ∧ (y˙ = −x)), but forbidding e.g. (x˙ = sin(y)) → m1 to avoid subtleties in the semantics of the formula. Additionally, Φ and the variables therein have the structure Φ = decl[0] ∧ · · · ∧ decl[k] ∧ init[0] ∧ trans[0, 1] ∧ . . . ∧ trans[(k − 1), k] ∧ target[k] (6) arising from the k-fold unwinding of the transition system, where decl[i] is the i-th instantiation of the system variables’ domain bounds; init[0] is the predicative encoding of the initial state applied to the 0-th variable instance, i.e. to the beginning of the trace; trans[i, i + 1] is the application of the transition predicate to the i-th and (i + 1)-th instances of the variables, e.g. instantiating a0 = a+1 to a[3] = a[2] + 1; and target[k] is the application of the target predicate to the last variable instance. ODE constraints occur only within the transition relation since they constrain the continuous flow behavior of the system. Satisfiability. As usual, a valuation σ, which maps each variable to a point from its domain, satisfies Φ iff the constraints satisfied under σ satisfy the Boolean structure of Φ. Satisfiability is straightforward for simple bounds and arithmetic constraints, but requires some explanation in the case of ODE constraints. ODEs occur only in the transition relation, which constrains the pre-post relation between any two successive instances of variables in a trace. Semantically, a trace is a sequence of snapshots of a real-time trajectory of the hybrid system. Hence, ODE constraints describe the behavior of the system between two such snapshots, i.e. describe trajectories emerging from the pre-valuation, following the dynamics described by the ODE, and finally reaching the post-valuation. A valuation σ thus satisfies a definitionally closed system of ODE constraints (each occurring variable itself being defined by one of the component ODE constraints), iff there exists a solution trajectory starting with the prevaluation and ending with the post-valuation, after a duration equal to the temporal length of the flow, as provided by the value of a special variable delta time. Example. Fig. 1 shows an example for a predicatively encoded hybrid system. In this example, the goal is to find a trajectory of the system that reaches a partially given system state. The system consists of two modes, which constrain the evolution of the continuous variable x by ODEs. Between these two modes, the system can switch arbitrarily at any point of time. Applying iSATODE to the encoding yields an unsatisfiability result for k = 0 and k = 1. That is, the solver proves that no trajectory with less than two transition steps reaches the target from the given initial state. Conversely, at least one mode switch is required. For two unwindings of the transition relation (k = 2), the solver finds a candidate solution box (Fig. 2). For each variable instance, this box contains a small interval for which the solver cannot prove
the absence of a valuation that satisfies the formula. While the algorithmic approach guarantees absence of solutions in case unsatisfiability is reported, candidate solutions may be spurious, i.e. they may or may not contain actual pointvalued solutions. If a box contains a solution, a k-step trace, it represents an actual trajectory over real time that passes through the steps of this trace. DECL int [0,1] mode; −− Discrete mode of the system. float [−100, 100] x; −− One continuous variable. float [0, 100] time; −− Global time and step duration float [0, 100] delta time; −− special variables. INIT x = 1.0; −− Initially, x = 1. time = 0; −− Start with time = 0. mode = 0; −− Start in mode = 0. TRANS mode = 0 −> (d.x / d.time = 2∗x); −− Dynamics described by mode = 1 −> (d.x / d.time = 0.2∗x); −− ODEs for each mode. −− Time progresses by time’ = time + delta time; delta time = 0 −> x’ = x; −− duration of step. TARGET x = 3.2; −− Want to reach x = 3.2 but no time >= 3.0; −− earlier than after time = 3.0.
Figure 1. Example encoding of a hybrid system. step 0 1 2
mode 0 1 –
delta time [0.051, 0.055] [5.281, 5.292] –
time [0, 0] [0.051, 0.055] [5.335, 5.343]
x [1, 1] [1.110, 1.114] [3.199, 3.201]
Figure 2. Candidate solution box. 3.2 Application to Set-Membership Estimation The predicative encoding of hybrid systems is a flexible tool, since it allows the integration of additional constraints into the formula, and thereby, to manipulate the space of solution trajectories. Such additional constraints may stem from measurements and error margins with which any solution trajectory shall be compatible. Together with the algorithmic guarantees that an unsatisfiability result is only generated, if no satisfying valuation exists, iSAT-ODE can be used to identify and prune off ranges of infeasible assignments to system parameters and state vectors. In practice, this means that iSAT will narrow down the domains for unknown system parameters and for the state vector based on the measurements supplied, using logical deduction to achieve such narrowing. The result either is refined knowledge about the parameter or state vector values that can explain the measurements, or—in case narrowing deduces empty intervals—a proof that the model is inconsistent with the measurements. Technically, this is achieved by first using only the deduction mechanisms of iSAT-ODE without any decisions by interval splitting to remove parts of the variables’ domains that cannot contain any solutions, hence by excluding those subranges for which iSAT-ODE has proven unsatisfiability of the formula. Second, a candidate solution box is searched with reactivated switching. In the third phase, the solver performs a bisection search for each of the variable bounds for which a range estimate is requested: while the distance between the bound of the remaining variable
442
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
range and the outer bound of the observed candidate solutions is still larger than the specified threshold, the variable is restricted to values between these two bounds, and a candidate solution is searched. If no such solution is found, the entire range can be excluded, since it is known not to contain any solution. Otherwise, the enclosure of the encountered candidate solutions is extended to cover this new box as well. 4. EXPERIMENTAL EVALUATION We demonstrate the applicability of our approach to different types of set membership estimation problems by the following experiments. 4.1 Parameter Estimation in a Continuous System First, we apply iSAT-ODE to a model of a chemical reaction, for which the basic kinematic dependencies between the involved variables are known, while the possible range of a reaction parameter is to be determined from an experimental measurement of the total reaction time. We model the reaction 2HI −→ H2 + I2 , and the assumption that the concentration of I2 and H2 depends quadratically on the concentration of HI by the following system of ODEs: c˙I2 = k · c2HI , c˙H2 = k · c2HI , c˙HI = −2 · k · c2HI . The parameter k is assumed to be constant, but unknown, and only coarsely enclosed by the interval [0, 10]. 1 As a first problem instance, we assume known initial concentrations cHI (0) = 1, cH2 (0) = cI2 (0) = 0 and an exact measurement of the time t = 5.7, when the concentration cHI (t) reaches 0.1. By adding this measurement as a target predicate, we obtain a model that is satisfiable by all valuations, for which the associated trajectories adhere to this measurement. Solving this system with iSAT-ODE, we obtain the following safe range estimate for the parameter k ∈ [0.77514648, 0.80000878]. Figure 3 shows two subsequently simulated trajectories for values from this range and intuitively supports this range estimate. 1
0.1035
simulated HI trajectory with k=0.776 simulated HI trajectory with k=0.800
0.9
0.103
0.8
0.1025
0.7
0.102
0.6
0.1015
0.5
0.101
0.4
0.1005
0.3
0.1
0.2
0.0995
0.1
As a second step, we extend this model to contain uncertainty in the initial condition and the measurement. We no longer assume the initial concentrations to be known exactly, but instead constrain cH2 (0), cI2 (0) ∈ [0, 0.02] and cHI (0) ∈ [0.95, 1]. The measurement is also uncertain; it is now assumed to have happened for t ∈ [5.6, 5.8], and the final concentration is taken as cHI (t) ∈ [0.08, 0.12]. The tightest range estimate we were able to generate with iSAT-ODE is given in Fig. 4. The most direct interpretation of these values is in terms of the variable instances of the formula comparable to the interpretation of a candidate solution box in Fig. 2 as a trace. However, the important difference is that here, the lines show the possible range of any solution, i.e. the system does not have a trajectory that does not pass through these two boxes and still satisfies the specified constraints. Positively, this means that—no matter which actual trajectory occurs— this range estimate guarantees that it passes through these boxes. Again, we can also extract a range estimate for the parameter k ∈ [0.622, 1.038] ∩ [0.612, 1.041]. The slightly different ranges for unwinding step 0 and 1 are caused by internally representing constant k in steps 0 and 1 by two different variables of the constraint formula, which is an artifact of the formula preprocessing of iSAT. Since we know that k is constant, we can safely intersect these two ranges. For the system without uncertainty, we have observed a runtime of 32 minutes on an Intel Core i7, 2.6 GHz processor for the range estimate of k, with a minimum splitting width (msw) down to which range refinements and interval splits in iSAT-ODE are performed of 0.01, and an absolute bound progress threshold (pr abs) of 0.001, which means that during the search for candidate solutions, new bounds are discarded, if they are closer than pr abs to already known bounds for the same variable. The solver runtimes for the range estimate reported in Fig. 4 for the system with uncertainty were significantly higher (56 hours), when using the same settings (except for disabling backwards propagation in the ODE-solver, which caused so significantly increased runtimes on the model that we aborted the solver run). However, increasing the msw and pr abs, hence accepting coarser results, allowed range estimates to be generated significantly faster. With an msw of 0.1 and pr abs of 0.01, we can enclose the same variables as in Fig. 4 within 8 hours, yielding e.g. k ∈ [0.505, 1.059].
0.099
0
0.0985
0
1
2
3
4
5
6
5.6
5.62
5.64
5.66
5.68
4.2 Set Membership Estimation for Hybrid Systems
5.7
Figure 3. Trajectories for the HI concentration over time with two values for the parameter k from the safe range estimate. (Obtained by numerical simulation.)
For x2 > k3 (mode above): √ x˙ 1 k1 − k2 x1 − x2 + k3 √ = √ x˙ 2 k2 x1 − x2 + k3 − k4 x2
k1 x1
step 0 1
k (0.622, 1.038) (0.612, 1.041)
cHI [0.949, 1.000] [0.079, 0.120]
cH2 [0.000, 0.021] [0.411, 0.489]
cI2 [0.000, 0.021] [0.409, 0.487]
·k2 k3
Figure 4. Range estimate for chemical reaction model with uncertainty.
x2
For x2 ≤ k3 (mode below): √ x˙ 1 k1 − k2 x1 = √ √ ·k4 x˙ 2 k2 x1 − k4 x2
Figure 5. Two tank system.
1
Since our focus is merely on the abstract problem structure, we do not claim chemical accuracy, and assume that an actual reaction parameter may well depend on additional factors like the temperature, which itself would be influenced by the reaction.
The second model is based on the two tank system from (Stursberg et al., 1997) used as a benchmark for hybrid system analysis tools. This system comprises two tanks
443
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
connected by a tube. The first tank has an inflow of constant k1 = 0.75 volume units, and its base is k3 = 0.5 length units above the base of the second tank. The connecting tube is characterized by a constant factor k2 = 1, which also characterizes the outflow of the system as k4 = 1. Figure 5 illustrates this system and formalizes the dynamic behavior of the liquid’s height x1 and x2 in the two tanks. The behavior switches between two dynamics, when x2 reaches the outlet from tank 1, and therefore, exerts a counter pressure against the incoming flow. Note that the model is bounded to the case x2 ≤ x1 + k3 , since it does not provide the dynamics for the inverse direction.
time
40 30 20 10 0
x1
0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25
x2
0.62 0.6 0.58 0.56 0.54 0.52 0.5 0.48
0
1
2
3
4
5
6
7
BMC depths
0
1
2
3
4
5
6
7
BMC depths
0
1
2
3
4
5
6
7
BMC depths
0
1
2
3
4
5
6
7
BMC depths
mode_below mode_above mode_dont_care
Remark: We also have to extend the model with a third mode, called mode dont care, which is located around a point on the switching surface where x˙ 2 = 0. For technical reasons, our encoding to enforce time progress is too weak at this point to avoid that the trajectory can immediately again perform a switch (hence chatter back and forth between the modes in this point infinitely often). By adding non-determinism around this point, but excluding the point itself from the reachable states, we still capture all possible system behavior—though with slightly decreased accuracy in case this small region is visited and the nondeterminism is used by the trajectory to jump contrary to the actual system dynamics—but gain the ability to find an upper bound of the number of required transition unwindings. Measuring at discrete time instants. Our first experiment with this system is based on three measurements of the variables x1 and x2 , at time t = 0.16, t = 0.48, and t = 5.28, which are assumed uncertain, e.g. x1 (0.16) ∈ [0.341− ε, 0.341 + ε], x2 (0.16) ∈ [0.557 − ε, 0.557 + ε]. The data were obtained from a simulated trajectory emerging from x1 (0) = 0.3, x2 (0) = 0.6. We modify the transition relation of the model in such a way that these measurements are interleaved with the normal transitions of the system, i.e., a transition can not only be a continuous flow following the dynamics of the current mode for some duration or a discrete mode switch, but can also be the “consumption” of a measurement event. Additionally, we enforce that a transition must have its maximum possible duration. That is, each step must end with reaching a predefined final point of time after the last measurement, or with reaching the switching surface from where a jump into the other mode becomes possible, or with reaching the point of time for the next measurement event. This modification is only aimed at ruling out solution traces that contain additional interruptions without a reason. The first goal is to calculate a number of unwindings of the transition relation, such that all measurements have been consumed, and all possible mode switches have been performed. We can now easily obtain this information by adding a target predicate that is true, when one of the measurements has not yet been consumed, or the final time is not yet reached. If this formula is satisfiable for a given unwinding depth, we know that such a trajectory still exists and increase the number of transition steps until the formula becomes unsatisfiable. For our model with the three measurements, iSAT-ODE can first prove unsatisfiability for depth 7. Removing the target predicate and instead adding the target that all measurements have been consumed, the model can be used to tighten the bounds of the initial values
Figure 6. Range estimate for the 2-tank system with seven unwindings (BMC depths) of the transition relation. x1 (0), x2 (0), which are only constrained to lie within the region where the system dynamics is defined. We report the range estimate for the variable instances of x1 , x2 , and the time variable in Fig. 6, along with the range estimates for the discrete mode in the bottom chart of the figure. Note that for step 3, the reconstructed mode is ambiguous: uncertainty makes feasible both modes mode below and mode above. Furthermore, like the structurally similar result in Fig. 4, this result must be interpreted in the following way: each trajectory of the system that is compatible with the measurements must pass through each of the boxes for steps 0 to 7. Especially, all trajectories must start in the initial box which has been determined to be x1 (0) ∈ [0.285, 0.315], x2 (0) ∈ [0.578, 0.623] and consequently start in mode above. This range estimate safely encloses the starting point of the simulated trajectory from which the measurements were obtained. Measuring height threshold. In our final experiment, instead of considering arbitrarily chosen measurements, we assume to have one sensor in each tank that measures when a certain height is crossed. Following the idea presented by Koutsoukos (2003) for a similar model of two tanks, placing one of these sensors at the level of switching, i.e. in our case observing when x2 (t) crosses k3 , allows direct observation of the time of mode changes. We can simplify our approach significantly, since most of the complexity in the previous experiment was caused by the necessity to determine a-priori unknown number of possible mode switches, and the consequently required number of transition steps. Starting from a basic encoding of the two tank system, we add a designated counter variable that is incremented in each step of the transition relation. Using this counter, we can iterate through the list of measurements and mode switches, directly enforcing constraints on the time to be reached in this step. While we assume the sensor to give an exact signal of whether the height level has been crossed, we allow for a small temporal uncertainty, i.e. do not enforce measurements to happen at a point of time, but instead within an interval around the assumed time of the observed level crossing. Figure 7 exemplifies the encoding of such an event and also shows that the a-priori knowledge of when events happen gives us the freedom to add additional interruptions of the flow at points of time, for which we know that no event occurs. When performing the range estimation, these additional steps allow better observation of the system state’s evolution. We illustrate
444
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
REFERENCES
step = 2 −> ( time’ >= 0.30 − EPS T −− EPS T is 0.01 and time’ <= 0.30 + EPS T and x2’ = 0.5 and (mode below’ <−> mode above)); step = 3 −> time’ = 0.40 ;
Figure 7. Part of the extended transition relation: encoding of a level crossing detection event leading to a mode switch and additional observation step. x1
0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25
x2
0.6
0
1
2
3
4
5
6
0
1
2
3
4
5
6
0
1
2
3
4
5
6
time
0.55 0.5 0.45 0.4 time
mode_below
mode_above
time
Figure 8. Range estimate for the measurements at discrete levels depicted over the valuation of the time variable. the range estimate in Figure 8, where we use the valuation of the time variable as x-axis. Each of the boxes drawn in this graph must be crossed by each trajectory satisfying all constraints. By having additional observation steps, where the time variable is forced to particular points, we now obtain actual enclosures: from the range estimate, we know that at time t = 0.4, the values for x1 and x2 must lie within the corresponding range estimate. Finally, we have also reconstructed the switching sequence. Runtimes. Checking unsatisfiability up to unwinding depth 7 in the first experiment has taken iSAT-ODE approximately 86 minutes on an AMD Opteron 8378, 2.4 GHz processor running 64 bit Linux using an msw of 0.01 and pr abs of 0.001. The range estimate for the 48 variable instances (six variables in eight instances) was completed within 5.7 hours using the same settings. The results reported in Figure 8, were obtained using an msw of 0.0625 and pr abs of 0.001 within 7.3 hours on the same machine. Using 0.1 as a coarser setting for the msw did not improve the runtime, whereas using a minimum splitting width of 0.01, we stopped the the solver after 27 hours. 5. CONCLUSION In this paper, we have shown that set membership estimation for both the discrete mode and the continuous state variables of uncertain and nonlinear hybrid dynamical systems can naturally be expressed as a SAT Modulo ODE formula. The latter are then efficiently solved using the satisfiability checker iSAT-ODE. The outcome of the experimental evaluations emphasizes the effectiveness of our approach. We expect the ongoing improvements of the underlying algorithms to also allow the application to higher dimensional hybrid systems in the future, hence improving the scalability of the approach.
Alur, R., Courcoubetis, C., Halbwachs, N., Henzinger, T., Ho, P.H., Nicollin, X., Olivero, A., Sifakis, J., and Yovine, S. (1995). The algorithmic analysis of hybrid systems. Theoretical Computer Science, 138, 3–34. Bemporad, A., Garulli, A., Paoletti, S., and Vicino, A. (2005). A bounded-error approach to piecewise affine system identification. Automatic Control, IEEE Transactions on, 50(10), 1567 – 1580. Benazera, E. and Trave-Massuyes, L. (2009). Set-theoretic estimation of hybrid system configurations. Systems, Man, and Cybernetics, Part B: Cybernetics, IEEE Transactions on, 39(5), 1277 –1291. Eggers, A., Ramdani, N., Nedialkov, N., and Fr¨ anzle, M. (2011). Improving SAT Modulo ODE for hybrid systems analysis by combining different enclosure methods. In G. Barthe, A. Pardo, and G. Schneider (eds.), Software Engineering and Formal Methods, volume 7041 of Lecture Notes in Computer Science, 172–187. Springer. Fr¨anzle, M., Herde, C., Ratschan, S., Schubert, T., and Teige, T. (2007). Efficient solving of large non-linear arithmetic constraint systems with complex boolean structure. JSAT Special Issue on Constraint Programming and SAT, 1, 209–236. Goldsztejn, A., Mullier, O., Eveillard, D., and Hosobe, H. (2010). Including ordinary differential equations based constraints in the standard CP framework. In D. Cohen (ed.), CP 2010, volume 6308 of Lecture Notes in Computer Science, 221–235. Springer. Jaulin, L. (2011). Range-only slam with occupancy maps: A set-membership approach. Robotics, IEEE Transactions on, 27(5), 1004 –1010. Kieffer, M. and Walter, E. (2011). Guaranteed estimation of the parameters of nonlinear continuous-time models: Contributions of interval analysis. International Journal of Adaptative Control and Signal Processing, 25, 191– 207. Koutsoukos, X. (2003). Estimation of hybrid systems using discrete sensors. In Decision and Control, 2003. Proc. 42nd IEEE Conf on, volume 1, 155–160. Le Bars, F., Sliwka, J., Jaulin, L., and Reynet, O. (2012). Set-membership state estimation with fleeting data. Automatica, 48(2), 381 – 387. Meslem, N. and Ramdani, N. (2011). Interval observer design based on nonlinear hybridization and practical stability analysis. International Journal of Adaptive Control and Signal Processing, 25(3), 228–248. Meslem, N., Ramdani, N., and Candau, Y. (2010). Using hybrid automata for set-membership state estimation with uncertain nonlinear continuous-time systems. Journal of Process Control, 4, 481–489. Milanese, M. and Novara, C. (2011). Unified set membership theory for identification, prediction and filtering of nonlinear systems. Automatica, 47(10), 2141 – 2151. Moisan, M., Bernard, O., and Gouz´e, J. (2009). Near optimal interval observers bundle for uncertain bioreactors. Automatica, 45, 291–295. Stursberg, O., Kowalewski, S., Hoffmann, I., and Preußig, J. (1997). Comparing timed and hybrid automata as approximations of continuous systems. In P. Antsaklis, W. Kohn, A. Nerode, and S. Sastry (eds.), Hybrid Systems IV, volume 1273 of LNCS, 361–377. Springer.
445
Set-Membership Estimation of Hybrid Systems via SAT Modulo ODE ? Andreas Eggers ∗ Nacim Ramdani ∗∗ Nedialko S. Nedialkov ∗∗∗ Martin Fr¨ anzle ∗ ∗
Carl von Ossietzky Universit¨ at, Oldenburg, Germany, (e-mail: {eggers—fraenzle}@informatik.uni-oldenburg.de). ∗∗ Universit´e d’Orl´eans, PRISME, 18020 Bourges, France, (e-mail: [email protected]) ∗∗∗ McMaster University, Hamilton, Ontario, Canada, (e-mail: [email protected]) Abstract: Set membership estimation (SME) of nonlinear hybrid systems is still a challenging issue. Although SME of nonlinear continuous systems has made significant progress recently, the direct extension of these methods to the hybrid case is not easy. Meanwhile, satisfiability (SAT) checkers for Boolean combinations of arithmetic constraints over real- and integer-valued variables have made significant progress, as they can effectively deal with algebraic constraints between variables and non-linear ODEs, what is denoted as SAT Modulo ODE. Finally, the corresponding solvers solve in a natural way the hybrid differential and algebraic constraints satisfaction problems that underlie SME of hybrid systems. This paper presents the application of such a SAT Modulo ODE solver to SME of hybrid dynamical systems. Keywords: Bounded error, constraint satisfaction problems, nonlinear systems, uncertain dynamic systems, state estimation, robust estimation, satisfiability. 1. INTRODUCTION State estimation is an important problem when addressing control or diagnosis issues of dynamical systems. Most controlled systems exhibit both smooth continuous dynamics and abrupt switches, hence can be efficiently modeled using hybrid automata, which combine discrete and continuous variables. In practice, the whole hybrid state remains only partially observable, as some switches and some continuous components cannot be directly observed. Hybrid state estimation aims at reconstructing both the discrete mode, hence the switching sequence, and the corresponding continuous state variables, based on a set of possibly discrete measurements, the knowledge of the hybrid model, and assumptions about the uncertainties and perturbations acting on the system. In many cases, it is more natural to assume that all uncertain quantities—measurement noise, model uncertainty and modelling errors—belong to a known set. This assumption is rather natural and requires much less data than statistical assumptions. In such an unknown-but-bounded-error (UBBE) framework, the estimation problem no longer has a unique solution, but there exists a set of state vectors that are consistent with measured data, the model structure and the prior error bounds. Set-membership estimation (SME) techniques allow the characterization of this solution set. SME tech? This work has been supported by the German Research Council DFG within SFB/TR 14 “Automatic Verification and Analysis of Complex Systems” (www.avacs.org), by the French National Research Agency under contract ANR 2011 INS 006 04 “MAGIC-SPS” (projects.laas.fr/ANR-MAGIC-SPS), and by the Natural Sciences and Engineering Research Council of Canada.
978-3-902823-06-9/12/$20.00 © 2012 IFAC
niques have reached a mature stage when dealing with linear or nonlinear discrete-time models; see a.o. (Milanese and Novara, 2011; Kieffer and Walter, 2011; Jaulin, 2011; Le Bars et al., 2012), and the references therein. Bemporad et al. (2005) addressed SME with piecewise affine systems. More recently, SME with continuous-time nonlinear systems has made significant progress; a.o. Moisan et al. (2009); Meslem et al. (2010); Meslem and Ramdani (2011) introduced new enclosing methods for nonlinear ordinary differential equations (ODE) for SME, and Goldsztejn et al. (2010) showed how to embed ODE into the numerical constraint satisfaction problem framework. SME with non linear hybrid systems is a challenging issue. Yet, nonlinear hybrid systems identification using point techniques is prohibitively time consuming, as it usually requires solving combinatorial nonlinear programming or optimization problems. The issue becomes even more difficult in set-membership framework, as SME solving techniques usually need to compute solution space coverings by means of multi-dimensional boxes. To the best of our knowledge, the only work addressing SME with nonlinear hybrid systems is by (Benazera and Trave-Massuyes, 2009), who addressed hybrid systems with discrete-time only nonlinear continuous dynamics. In this paper, we address SME of nonlinear hybrid systems, whose continuous dynamics are described using nonlinear ODEs. We show that SME of hybrid systems boils down to solving differential and algebraic constraints satisfaction problems, which in turn can be written as SAT Modulo ODE formulas, hence solved using appropriate satisfiability checkers— notions to be defined in following sections.
440
10.3182/20120711-3-BE-2027.00292
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
Outline. Sect. 2 recalls hybrid system modeling and SME. Sect. 3 introduces the main contribution of this work, namely the use of satisfiability checkers for solving SME with hybrid systems. Sect. 4 gathers experimental evaluation on two examples, including a hybrid system. 2. SET-MEMBERSHIP ESTIMATION OF HYBRID SYSTEMS We consider a hybrid automaton (Alur et al., 1995) given by H = (Q, D, P, Σ, A, Inv, F). (1) Here Q is a set of locations. Given a location q ∈ Q, the continuous dynamics, and hence flow transitions, are described by non-autonomous differential equations fq ∈ F of the form flow(q) : x(t) ˙ = fq (x, p, t), (2) where fq : D × P × R+ 7→ D is nonlinear and assumed sufficiently smooth over D ⊆ Rn , with dimension n that may depend on q, and p ∈ P, where P is an uncertainty domain for the parameter vector p. Inv is an invariant, which assigns a domain to the continuous state space of each location. It is defined by the following system of inequalities: Inv(q) : νq (x(t), p, t) < 0, (3) where νq : D × P × R+ 7→ Rm is a vector-valued nonlinear function, the negativity constraint applies componentwise, and the number m of inequalities may also depend on q. The set A is the set of discrete transitions {e = (q → q 0 )} each given by the 5-uple (q, guard, σ, ρ, q 0 ), where q and q 0 represent upstream and downstream locations respectively; guard is a condition given as the system of equations guard(e) : γe (x(t), p, t) = 0, (4) σ is an event, and ρ is a reset function. A transition q → q 0 may occur when the continuous flow resides within the guard set, i.e. when the continuous state satisfies condition (4). Let us also consider the following measurement equation y(t) = gq (x(t), p, t), (5) where function gq : D × P × R+ 7→ Rm is nonlinear. We assume that one has access to measurements of the output vector y(t). The latter is usually measured at some sampling time instants tj taken over a time-grid t0 < t1 < t2 < . . . < tnT . The sampling intervals need not be constant, as it may be the case when measuring for instance reactant concentrations in biochemical reactors. The output vector may also be measured using discrete sensors, i.e. sensors that monitor only a single threshold value such as critical fluid level in chemical processes or proximity sensors in robotics applications. Discrete sensors are suitable for monitoring guard conditions that trigger discrete transitions in hybrid systems as emphasized by Koutsoukos (2003). Measurement noise is taken as a discrete time signal assumed additive and bounded with known bounds. Denote by y i the feasible domain for output vector yi at time ti . Recall that the parameter vector p is uncertain and lay within a bounded set with known bounds P. We assume that the initial hybrid state vector is only partially known, i.e. the actual discrete mode may be unknown or the continuous component taken in a bounded domain x(t0 ) ∈ x0 .
The goal of the bounded error estimation is to compute conservative outer enclosures for feasible sets for both the discrete modes and associated continuous variables that are consistent with the feasible domains for measurements, the chosen hybrid model, and the assumptions about the uncertainties and perturbations acting on the system. That is, given two measurements y i and y i+1 gathered at the two time instants ti and ti+1 , the estimation problem boils down to simultaneously: (i) Detect and identify all possible discrete transitions (t? , e) that may occur at time t? ∈ [ti , ti+1 ]. Because there are uncertainties in either initial state vector, continuous dynamics or measurements, there might be a continuum of time instants, where events may occur. There may also exist several such time intervals in [ti , ti+1 ]. (ii) Reconstruct the sets of hybrid states (q(ti ), x(ti )) and (q(ti+1 ), x(ti+1 )) that are consistent with the switching sequence reconstructed in (i), the modelling (1)-(5), and the inclusions y(ti ) ∈ y i and y(ti+1 ) ∈ y i+1 . These two steps must be done for all i, 0 ≤ i ≤ nT − 1. Finally, this set membership estimation problem for hybrid systems can naturally be stated as finding feasible solution sets for Boolean combinations of both arithmetic and differential constraints over real- and integer-valued variables. This is the purpose of the iSAT-ODE solver we introduce in the next section. 3. THE ISAT-ODE SOLVER In (Eggers et al., 2011), we presented an extension of the iSAT solver (Fr¨anzle et al., 2007) with a combination of enclosure methods for ODEs. In this section we concentrate on the application of iSAT-ODE on the set estimation problem and refer the reader to the above references for algorithmic details. 3.1 Encoding and Analysis of Hybrid Systems Our approach to constraint-based analysis of hybrid systems is based on a predicative encoding of the hybrid system’s trajectories. This encoding consists of a predicate describing the valid initial states, a k-fold unwinding of a formula connecting pairs of states by continuous and discrete transitions and—when considering classical bounded model checking problems—a target predicate that describes the final state of the system, whose reachability is to be analyzed. SAT Modulo ODE. Let Φ be a SAT Modulo ODE formula. This is a quantifier-free Boolean combination of arithmetic constraints over bounded real-, integer-, and Booleanvalued variables, simple bounds, and ODE constraints over real variables with the following properties:
441
• arithmetic constraints over variables x, y, and z are of the form x ∼ ◦(y, z) or x ∼ ◦(y), where ∼ is a relational operator from {<, ≤, =, ≥, >}, and ◦ is a total unary or binary operator from {+, −, ·, sin, cos, powN , exp, min, max}; • simple bounds are of the form x ∼ c with ∼ as above a relational operator, x a variable, and c a constant; and
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
• ODE constraints are sufficiently smooth, time invariant, and given by x˙ i = dxi /dt = f (x1 , . . . , xn ) with all occurring variables xi themselves being defined by ODE constraints and f being a function com√ posed of {+, −, ·, /, powN , exp, ln, , sin, cos}. These ODE constraints must occur only under an even number of negations in the formula, allowing e.g. m1 → ((x˙ = sin(y)) ∧ (y˙ = −x)), but forbidding e.g. (x˙ = sin(y)) → m1 to avoid subtleties in the semantics of the formula. Additionally, Φ and the variables therein have the structure Φ = decl[0] ∧ · · · ∧ decl[k] ∧ init[0] ∧ trans[0, 1] ∧ . . . ∧ trans[(k − 1), k] ∧ target[k] (6) arising from the k-fold unwinding of the transition system, where decl[i] is the i-th instantiation of the system variables’ domain bounds; init[0] is the predicative encoding of the initial state applied to the 0-th variable instance, i.e. to the beginning of the trace; trans[i, i + 1] is the application of the transition predicate to the i-th and (i + 1)-th instances of the variables, e.g. instantiating a0 = a+1 to a[3] = a[2] + 1; and target[k] is the application of the target predicate to the last variable instance. ODE constraints occur only within the transition relation since they constrain the continuous flow behavior of the system. Satisfiability. As usual, a valuation σ, which maps each variable to a point from its domain, satisfies Φ iff the constraints satisfied under σ satisfy the Boolean structure of Φ. Satisfiability is straightforward for simple bounds and arithmetic constraints, but requires some explanation in the case of ODE constraints. ODEs occur only in the transition relation, which constrains the pre-post relation between any two successive instances of variables in a trace. Semantically, a trace is a sequence of snapshots of a real-time trajectory of the hybrid system. Hence, ODE constraints describe the behavior of the system between two such snapshots, i.e. describe trajectories emerging from the pre-valuation, following the dynamics described by the ODE, and finally reaching the post-valuation. A valuation σ thus satisfies a definitionally closed system of ODE constraints (each occurring variable itself being defined by one of the component ODE constraints), iff there exists a solution trajectory starting with the prevaluation and ending with the post-valuation, after a duration equal to the temporal length of the flow, as provided by the value of a special variable delta time. Example. Fig. 1 shows an example for a predicatively encoded hybrid system. In this example, the goal is to find a trajectory of the system that reaches a partially given system state. The system consists of two modes, which constrain the evolution of the continuous variable x by ODEs. Between these two modes, the system can switch arbitrarily at any point of time. Applying iSATODE to the encoding yields an unsatisfiability result for k = 0 and k = 1. That is, the solver proves that no trajectory with less than two transition steps reaches the target from the given initial state. Conversely, at least one mode switch is required. For two unwindings of the transition relation (k = 2), the solver finds a candidate solution box (Fig. 2). For each variable instance, this box contains a small interval for which the solver cannot prove
the absence of a valuation that satisfies the formula. While the algorithmic approach guarantees absence of solutions in case unsatisfiability is reported, candidate solutions may be spurious, i.e. they may or may not contain actual pointvalued solutions. If a box contains a solution, a k-step trace, it represents an actual trajectory over real time that passes through the steps of this trace. DECL int [0,1] mode; −− Discrete mode of the system. float [−100, 100] x; −− One continuous variable. float [0, 100] time; −− Global time and step duration float [0, 100] delta time; −− special variables. INIT x = 1.0; −− Initially, x = 1. time = 0; −− Start with time = 0. mode = 0; −− Start in mode = 0. TRANS mode = 0 −> (d.x / d.time = 2∗x); −− Dynamics described by mode = 1 −> (d.x / d.time = 0.2∗x); −− ODEs for each mode. −− Time progresses by time’ = time + delta time; delta time = 0 −> x’ = x; −− duration of step. TARGET x = 3.2; −− Want to reach x = 3.2 but no time >= 3.0; −− earlier than after time = 3.0.
Figure 1. Example encoding of a hybrid system. step 0 1 2
mode 0 1 –
delta time [0.051, 0.055] [5.281, 5.292] –
time [0, 0] [0.051, 0.055] [5.335, 5.343]
x [1, 1] [1.110, 1.114] [3.199, 3.201]
Figure 2. Candidate solution box. 3.2 Application to Set-Membership Estimation The predicative encoding of hybrid systems is a flexible tool, since it allows the integration of additional constraints into the formula, and thereby, to manipulate the space of solution trajectories. Such additional constraints may stem from measurements and error margins with which any solution trajectory shall be compatible. Together with the algorithmic guarantees that an unsatisfiability result is only generated, if no satisfying valuation exists, iSAT-ODE can be used to identify and prune off ranges of infeasible assignments to system parameters and state vectors. In practice, this means that iSAT will narrow down the domains for unknown system parameters and for the state vector based on the measurements supplied, using logical deduction to achieve such narrowing. The result either is refined knowledge about the parameter or state vector values that can explain the measurements, or—in case narrowing deduces empty intervals—a proof that the model is inconsistent with the measurements. Technically, this is achieved by first using only the deduction mechanisms of iSAT-ODE without any decisions by interval splitting to remove parts of the variables’ domains that cannot contain any solutions, hence by excluding those subranges for which iSAT-ODE has proven unsatisfiability of the formula. Second, a candidate solution box is searched with reactivated switching. In the third phase, the solver performs a bisection search for each of the variable bounds for which a range estimate is requested: while the distance between the bound of the remaining variable
442
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
range and the outer bound of the observed candidate solutions is still larger than the specified threshold, the variable is restricted to values between these two bounds, and a candidate solution is searched. If no such solution is found, the entire range can be excluded, since it is known not to contain any solution. Otherwise, the enclosure of the encountered candidate solutions is extended to cover this new box as well. 4. EXPERIMENTAL EVALUATION We demonstrate the applicability of our approach to different types of set membership estimation problems by the following experiments. 4.1 Parameter Estimation in a Continuous System First, we apply iSAT-ODE to a model of a chemical reaction, for which the basic kinematic dependencies between the involved variables are known, while the possible range of a reaction parameter is to be determined from an experimental measurement of the total reaction time. We model the reaction 2HI −→ H2 + I2 , and the assumption that the concentration of I2 and H2 depends quadratically on the concentration of HI by the following system of ODEs: c˙I2 = k · c2HI , c˙H2 = k · c2HI , c˙HI = −2 · k · c2HI . The parameter k is assumed to be constant, but unknown, and only coarsely enclosed by the interval [0, 10]. 1 As a first problem instance, we assume known initial concentrations cHI (0) = 1, cH2 (0) = cI2 (0) = 0 and an exact measurement of the time t = 5.7, when the concentration cHI (t) reaches 0.1. By adding this measurement as a target predicate, we obtain a model that is satisfiable by all valuations, for which the associated trajectories adhere to this measurement. Solving this system with iSAT-ODE, we obtain the following safe range estimate for the parameter k ∈ [0.77514648, 0.80000878]. Figure 3 shows two subsequently simulated trajectories for values from this range and intuitively supports this range estimate. 1
0.1035
simulated HI trajectory with k=0.776 simulated HI trajectory with k=0.800
0.9
0.103
0.8
0.1025
0.7
0.102
0.6
0.1015
0.5
0.101
0.4
0.1005
0.3
0.1
0.2
0.0995
0.1
As a second step, we extend this model to contain uncertainty in the initial condition and the measurement. We no longer assume the initial concentrations to be known exactly, but instead constrain cH2 (0), cI2 (0) ∈ [0, 0.02] and cHI (0) ∈ [0.95, 1]. The measurement is also uncertain; it is now assumed to have happened for t ∈ [5.6, 5.8], and the final concentration is taken as cHI (t) ∈ [0.08, 0.12]. The tightest range estimate we were able to generate with iSAT-ODE is given in Fig. 4. The most direct interpretation of these values is in terms of the variable instances of the formula comparable to the interpretation of a candidate solution box in Fig. 2 as a trace. However, the important difference is that here, the lines show the possible range of any solution, i.e. the system does not have a trajectory that does not pass through these two boxes and still satisfies the specified constraints. Positively, this means that—no matter which actual trajectory occurs— this range estimate guarantees that it passes through these boxes. Again, we can also extract a range estimate for the parameter k ∈ [0.622, 1.038] ∩ [0.612, 1.041]. The slightly different ranges for unwinding step 0 and 1 are caused by internally representing constant k in steps 0 and 1 by two different variables of the constraint formula, which is an artifact of the formula preprocessing of iSAT. Since we know that k is constant, we can safely intersect these two ranges. For the system without uncertainty, we have observed a runtime of 32 minutes on an Intel Core i7, 2.6 GHz processor for the range estimate of k, with a minimum splitting width (msw) down to which range refinements and interval splits in iSAT-ODE are performed of 0.01, and an absolute bound progress threshold (pr abs) of 0.001, which means that during the search for candidate solutions, new bounds are discarded, if they are closer than pr abs to already known bounds for the same variable. The solver runtimes for the range estimate reported in Fig. 4 for the system with uncertainty were significantly higher (56 hours), when using the same settings (except for disabling backwards propagation in the ODE-solver, which caused so significantly increased runtimes on the model that we aborted the solver run). However, increasing the msw and pr abs, hence accepting coarser results, allowed range estimates to be generated significantly faster. With an msw of 0.1 and pr abs of 0.01, we can enclose the same variables as in Fig. 4 within 8 hours, yielding e.g. k ∈ [0.505, 1.059].
0.099
0
0.0985
0
1
2
3
4
5
6
5.6
5.62
5.64
5.66
5.68
4.2 Set Membership Estimation for Hybrid Systems
5.7
Figure 3. Trajectories for the HI concentration over time with two values for the parameter k from the safe range estimate. (Obtained by numerical simulation.)
For x2 > k3 (mode above): √ x˙ 1 k1 − k2 x1 − x2 + k3 √ = √ x˙ 2 k2 x1 − x2 + k3 − k4 x2
k1 x1
step 0 1
k (0.622, 1.038) (0.612, 1.041)
cHI [0.949, 1.000] [0.079, 0.120]
cH2 [0.000, 0.021] [0.411, 0.489]
cI2 [0.000, 0.021] [0.409, 0.487]
·k2 k3
Figure 4. Range estimate for chemical reaction model with uncertainty.
x2
For x2 ≤ k3 (mode below): √ x˙ 1 k1 − k2 x1 = √ √ ·k4 x˙ 2 k2 x1 − k4 x2
Figure 5. Two tank system.
1
Since our focus is merely on the abstract problem structure, we do not claim chemical accuracy, and assume that an actual reaction parameter may well depend on additional factors like the temperature, which itself would be influenced by the reaction.
The second model is based on the two tank system from (Stursberg et al., 1997) used as a benchmark for hybrid system analysis tools. This system comprises two tanks
443
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
connected by a tube. The first tank has an inflow of constant k1 = 0.75 volume units, and its base is k3 = 0.5 length units above the base of the second tank. The connecting tube is characterized by a constant factor k2 = 1, which also characterizes the outflow of the system as k4 = 1. Figure 5 illustrates this system and formalizes the dynamic behavior of the liquid’s height x1 and x2 in the two tanks. The behavior switches between two dynamics, when x2 reaches the outlet from tank 1, and therefore, exerts a counter pressure against the incoming flow. Note that the model is bounded to the case x2 ≤ x1 + k3 , since it does not provide the dynamics for the inverse direction.
time
40 30 20 10 0
x1
0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25
x2
0.62 0.6 0.58 0.56 0.54 0.52 0.5 0.48
0
1
2
3
4
5
6
7
BMC depths
0
1
2
3
4
5
6
7
BMC depths
0
1
2
3
4
5
6
7
BMC depths
0
1
2
3
4
5
6
7
BMC depths
mode_below mode_above mode_dont_care
Remark: We also have to extend the model with a third mode, called mode dont care, which is located around a point on the switching surface where x˙ 2 = 0. For technical reasons, our encoding to enforce time progress is too weak at this point to avoid that the trajectory can immediately again perform a switch (hence chatter back and forth between the modes in this point infinitely often). By adding non-determinism around this point, but excluding the point itself from the reachable states, we still capture all possible system behavior—though with slightly decreased accuracy in case this small region is visited and the nondeterminism is used by the trajectory to jump contrary to the actual system dynamics—but gain the ability to find an upper bound of the number of required transition unwindings. Measuring at discrete time instants. Our first experiment with this system is based on three measurements of the variables x1 and x2 , at time t = 0.16, t = 0.48, and t = 5.28, which are assumed uncertain, e.g. x1 (0.16) ∈ [0.341− ε, 0.341 + ε], x2 (0.16) ∈ [0.557 − ε, 0.557 + ε]. The data were obtained from a simulated trajectory emerging from x1 (0) = 0.3, x2 (0) = 0.6. We modify the transition relation of the model in such a way that these measurements are interleaved with the normal transitions of the system, i.e., a transition can not only be a continuous flow following the dynamics of the current mode for some duration or a discrete mode switch, but can also be the “consumption” of a measurement event. Additionally, we enforce that a transition must have its maximum possible duration. That is, each step must end with reaching a predefined final point of time after the last measurement, or with reaching the switching surface from where a jump into the other mode becomes possible, or with reaching the point of time for the next measurement event. This modification is only aimed at ruling out solution traces that contain additional interruptions without a reason. The first goal is to calculate a number of unwindings of the transition relation, such that all measurements have been consumed, and all possible mode switches have been performed. We can now easily obtain this information by adding a target predicate that is true, when one of the measurements has not yet been consumed, or the final time is not yet reached. If this formula is satisfiable for a given unwinding depth, we know that such a trajectory still exists and increase the number of transition steps until the formula becomes unsatisfiable. For our model with the three measurements, iSAT-ODE can first prove unsatisfiability for depth 7. Removing the target predicate and instead adding the target that all measurements have been consumed, the model can be used to tighten the bounds of the initial values
Figure 6. Range estimate for the 2-tank system with seven unwindings (BMC depths) of the transition relation. x1 (0), x2 (0), which are only constrained to lie within the region where the system dynamics is defined. We report the range estimate for the variable instances of x1 , x2 , and the time variable in Fig. 6, along with the range estimates for the discrete mode in the bottom chart of the figure. Note that for step 3, the reconstructed mode is ambiguous: uncertainty makes feasible both modes mode below and mode above. Furthermore, like the structurally similar result in Fig. 4, this result must be interpreted in the following way: each trajectory of the system that is compatible with the measurements must pass through each of the boxes for steps 0 to 7. Especially, all trajectories must start in the initial box which has been determined to be x1 (0) ∈ [0.285, 0.315], x2 (0) ∈ [0.578, 0.623] and consequently start in mode above. This range estimate safely encloses the starting point of the simulated trajectory from which the measurements were obtained. Measuring height threshold. In our final experiment, instead of considering arbitrarily chosen measurements, we assume to have one sensor in each tank that measures when a certain height is crossed. Following the idea presented by Koutsoukos (2003) for a similar model of two tanks, placing one of these sensors at the level of switching, i.e. in our case observing when x2 (t) crosses k3 , allows direct observation of the time of mode changes. We can simplify our approach significantly, since most of the complexity in the previous experiment was caused by the necessity to determine a-priori unknown number of possible mode switches, and the consequently required number of transition steps. Starting from a basic encoding of the two tank system, we add a designated counter variable that is incremented in each step of the transition relation. Using this counter, we can iterate through the list of measurements and mode switches, directly enforcing constraints on the time to be reached in this step. While we assume the sensor to give an exact signal of whether the height level has been crossed, we allow for a small temporal uncertainty, i.e. do not enforce measurements to happen at a point of time, but instead within an interval around the assumed time of the observed level crossing. Figure 7 exemplifies the encoding of such an event and also shows that the a-priori knowledge of when events happen gives us the freedom to add additional interruptions of the flow at points of time, for which we know that no event occurs. When performing the range estimation, these additional steps allow better observation of the system state’s evolution. We illustrate
444
16th IFAC Symposium on System Identification Brussels, Belgium. July 11-13, 2012
REFERENCES
step = 2 −> ( time’ >= 0.30 − EPS T −− EPS T is 0.01 and time’ <= 0.30 + EPS T and x2’ = 0.5 and (mode below’ <−> mode above)); step = 3 −> time’ = 0.40 ;
Figure 7. Part of the extended transition relation: encoding of a level crossing detection event leading to a mode switch and additional observation step. x1
0.65 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25
x2
0.6
0
1
2
3
4
5
6
0
1
2
3
4
5
6
0
1
2
3
4
5
6
time
0.55 0.5 0.45 0.4 time
mode_below
mode_above
time
Figure 8. Range estimate for the measurements at discrete levels depicted over the valuation of the time variable. the range estimate in Figure 8, where we use the valuation of the time variable as x-axis. Each of the boxes drawn in this graph must be crossed by each trajectory satisfying all constraints. By having additional observation steps, where the time variable is forced to particular points, we now obtain actual enclosures: from the range estimate, we know that at time t = 0.4, the values for x1 and x2 must lie within the corresponding range estimate. Finally, we have also reconstructed the switching sequence. Runtimes. Checking unsatisfiability up to unwinding depth 7 in the first experiment has taken iSAT-ODE approximately 86 minutes on an AMD Opteron 8378, 2.4 GHz processor running 64 bit Linux using an msw of 0.01 and pr abs of 0.001. The range estimate for the 48 variable instances (six variables in eight instances) was completed within 5.7 hours using the same settings. The results reported in Figure 8, were obtained using an msw of 0.0625 and pr abs of 0.001 within 7.3 hours on the same machine. Using 0.1 as a coarser setting for the msw did not improve the runtime, whereas using a minimum splitting width of 0.01, we stopped the the solver after 27 hours. 5. CONCLUSION In this paper, we have shown that set membership estimation for both the discrete mode and the continuous state variables of uncertain and nonlinear hybrid dynamical systems can naturally be expressed as a SAT Modulo ODE formula. The latter are then efficiently solved using the satisfiability checker iSAT-ODE. The outcome of the experimental evaluations emphasizes the effectiveness of our approach. We expect the ongoing improvements of the underlying algorithms to also allow the application to higher dimensional hybrid systems in the future, hence improving the scalability of the approach.
Alur, R., Courcoubetis, C., Halbwachs, N., Henzinger, T., Ho, P.H., Nicollin, X., Olivero, A., Sifakis, J., and Yovine, S. (1995). The algorithmic analysis of hybrid systems. Theoretical Computer Science, 138, 3–34. Bemporad, A., Garulli, A., Paoletti, S., and Vicino, A. (2005). A bounded-error approach to piecewise affine system identification. Automatic Control, IEEE Transactions on, 50(10), 1567 – 1580. Benazera, E. and Trave-Massuyes, L. (2009). Set-theoretic estimation of hybrid system configurations. Systems, Man, and Cybernetics, Part B: Cybernetics, IEEE Transactions on, 39(5), 1277 –1291. Eggers, A., Ramdani, N., Nedialkov, N., and Fr¨ anzle, M. (2011). Improving SAT Modulo ODE for hybrid systems analysis by combining different enclosure methods. In G. Barthe, A. Pardo, and G. Schneider (eds.), Software Engineering and Formal Methods, volume 7041 of Lecture Notes in Computer Science, 172–187. Springer. Fr¨anzle, M., Herde, C., Ratschan, S., Schubert, T., and Teige, T. (2007). Efficient solving of large non-linear arithmetic constraint systems with complex boolean structure. JSAT Special Issue on Constraint Programming and SAT, 1, 209–236. Goldsztejn, A., Mullier, O., Eveillard, D., and Hosobe, H. (2010). Including ordinary differential equations based constraints in the standard CP framework. In D. Cohen (ed.), CP 2010, volume 6308 of Lecture Notes in Computer Science, 221–235. Springer. Jaulin, L. (2011). Range-only slam with occupancy maps: A set-membership approach. Robotics, IEEE Transactions on, 27(5), 1004 –1010. Kieffer, M. and Walter, E. (2011). Guaranteed estimation of the parameters of nonlinear continuous-time models: Contributions of interval analysis. International Journal of Adaptative Control and Signal Processing, 25, 191– 207. Koutsoukos, X. (2003). Estimation of hybrid systems using discrete sensors. In Decision and Control, 2003. Proc. 42nd IEEE Conf on, volume 1, 155–160. Le Bars, F., Sliwka, J., Jaulin, L., and Reynet, O. (2012). Set-membership state estimation with fleeting data. Automatica, 48(2), 381 – 387. Meslem, N. and Ramdani, N. (2011). Interval observer design based on nonlinear hybridization and practical stability analysis. International Journal of Adaptive Control and Signal Processing, 25(3), 228–248. Meslem, N., Ramdani, N., and Candau, Y. (2010). Using hybrid automata for set-membership state estimation with uncertain nonlinear continuous-time systems. Journal of Process Control, 4, 481–489. Milanese, M. and Novara, C. (2011). Unified set membership theory for identification, prediction and filtering of nonlinear systems. Automatica, 47(10), 2141 – 2151. Moisan, M., Bernard, O., and Gouz´e, J. (2009). Near optimal interval observers bundle for uncertain bioreactors. Automatica, 45, 291–295. Stursberg, O., Kowalewski, S., Hoffmann, I., and Preußig, J. (1997). Comparing timed and hybrid automata as approximations of continuous systems. In P. Antsaklis, W. Kohn, A. Nerode, and S. Sastry (eds.), Hybrid Systems IV, volume 1273 of LNCS, 361–377. Springer.
445