Scalable computation of viability kernels and a viability-theoretic approach to guaranteeing safety for closed-loop medical devices by John Norman Maidens BSc. (Honors) in Mathematics, University of Alberta, 2010 a thesis submitted in partial fulfillment of the requirements for the degree of Master of Applied Science in the faculty of graduate studies (Biomedical Engineering) The University Of British Columbia (Vancouver) July 2012 c© John Norman Maidens, 2012 Abstract As closed-loop controllers become increasingly prevalent in medical tech- nology, increasing emphasis is being placed on ensuring that such systems operate in a safe manner. In our approach to guaranteeing the safe operation of a physiologic closed-loop control system, we wish to provide a mathemat- ical guarantee that, despite limited control authority, the system’s state can be confined to a region designated as safe. The largest subset of the safe region for which there exists an admissible control input that keeps the state within the safe region is known as the viability kernel, or maximal controlled invariant set. Many methods are known for computing viability kernels in low-dimensional systems, but these existing methods rely on gridding the state space and hence their time complexity increases exponentially with the state dimen- sion. In this thesis we describe a new connection between reachability and viability theory that enables us to approximate the viability kernel using Lagrangian methods which scale well with the state dimension. We present four new viability kernel approximation algorithms using polytope-, ellipsoid- and support vector-based set representations and we compare their performances in terms of accuracy and scalability with the state dimension. Using the support vector and ellipsoidal techniques, we are able to accurately approximate the viability kernel for systems of much larger state dimension than was previously feasible using existing Eulerian methods. We also present a viability theoretic solution to the problem of determin- ing when a physiologic closed-loop control system should initiate a fallback ii mode of operation. The viability-based method allows impending safety vio- lations to be detected in advance, allowing the fallback mode to be initiated earlier than using a naive approach. Our new approach to fallback mode initiation is examined in two sample contexts: the closed-loop control of carbon dioxide partial pressure under mechanical ventilation, and the con- trol of the concentration of the anaesthetic drug Propofol using a paediatric model of Propofol pharmacokinetics. iii Preface The research presented in this thesis is the result of a close and fruitful col- laboration with Shahab Kaynama. Early in my Master’s program, Shahab and I began to discuss the idea of using reachability techniques to approxi- mate viability kernels, something that Shahab had been thinking about for a few months, but had not yet been able to formulate in a satisfactory way. I was responsible for the first results for continuous-time systems given here in Theorems 4 and 5 and Shahab provided the extension to discrete-time systems presented here as Theorem 3. Shahab developed the first of the Lagrangian algorithms (Algorithms 3.1 and 3.3) while I was responsible for developing the approximate polytopic algorithm (Algorithm 3.2) and the support vector algorithm (Algorithm 3.4). The research described in this thesis has resulted in two publications. A paper entitled “Computing the viability kernel using maximal reachable sets” was drafted by Shahab with myself, Ian M. Mitchell, Meeko Oishi and Guy A. Dumont as co-authors and it was presented at Hybrid Systems: Computation and Control 2012 in Beijing. I have drafted a manuscript en- titled “Lagrangian methods for approximating the viability kernel in high- dimensional systems” with Shahab Kaynama, Ian M. Mitchell, Meeko Oishi and Guy A. Dumont as co-authors, which is currently under review. Mate- rial from this second publication forms much of the body of Chapters 2 and 3 of this thesis. The material appearing in Chapter 4 is more recent, and has not been submitted for publication. iv Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . v List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . xii 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Computational approaches to verification . . . . . . . . . . . 2 1.1.1 Simulation . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.2 Barrier functions and invariant sets . . . . . . . . . . . 3 1.1.3 Reachability-based approaches . . . . . . . . . . . . . 4 1.2 Viability kernel approximation methods . . . . . . . . . . . . 4 1.2.1 Viability kernel algorithm . . . . . . . . . . . . . . . . 5 1.2.2 Supervised classification . . . . . . . . . . . . . . . . . 5 1.2.3 Approximate dynamic programming . . . . . . . . . . 6 1.2.4 Simulated annealing . . . . . . . . . . . . . . . . . . . 6 1.3 Contributions of this thesis . . . . . . . . . . . . . . . . . . . 7 2 Establishing connections between viability and reachability 9 2.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.1.1 Elements of viability theory . . . . . . . . . . . . . . . 11 v 2.1.2 Reachability under input constraints . . . . . . . . . . 13 2.2 The Viability Kernel Algorithm . . . . . . . . . . . . . . . . . 14 2.3 Computing viability kernels using reachability techniques . . 17 2.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3 Lagrangian algorithms for approximating viability kernels in linear systems . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.1 Set representations for linear reachability . . . . . . . . . . . 26 3.1.1 Convex polytopes . . . . . . . . . . . . . . . . . . . . . 26 3.1.2 Ellipsoids . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.1.3 Support functions . . . . . . . . . . . . . . . . . . . . 28 3.1.4 Support vectors . . . . . . . . . . . . . . . . . . . . . . 29 3.2 Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.2.1 Exact polytopic method . . . . . . . . . . . . . . . . . 32 3.2.2 Approximate polytopic method . . . . . . . . . . . . . 33 3.2.3 Ellipsoidal method . . . . . . . . . . . . . . . . . . . . 35 3.2.4 Support vector method . . . . . . . . . . . . . . . . . 36 3.3 Comparison of algorithms . . . . . . . . . . . . . . . . . . . . 40 3.3.1 Accuracy . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.3.2 Scalability . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.4 High-dimensional example . . . . . . . . . . . . . . . . . . . . 49 3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4 Fallback mode initiation for physiologic closed-loop con- trollers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 4.1 Constraints for physiologic variables . . . . . . . . . . . . . . 54 4.1.1 Input and state constraints . . . . . . . . . . . . . . . 54 4.1.2 Refining constraints for early fallback mode initiation 54 4.2 Illustrative example: closed-loop mechanical ventilation . . . 55 4.3 Closed-loop control of anaesthesia . . . . . . . . . . . . . . . 60 4.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 5.1 Future research . . . . . . . . . . . . . . . . . . . . . . . . . . 63 vi 5.1.1 Differential games . . . . . . . . . . . . . . . . . . . . 63 5.1.2 Improved computations with support vectors . . . . . 64 5.1.3 Time-scale separation methods . . . . . . . . . . . . . 64 5.1.4 Application to anaesthesia . . . . . . . . . . . . . . . . 65 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 A Proof of Theorem 7 . . . . . . . . . . . . . . . . . . . . . . . . 70 A.1 Convex conjugation and duality . . . . . . . . . . . . . . . . . 70 A.2 Subgradients and support vectors . . . . . . . . . . . . . . . . 71 A.3 Theorem 7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 vii List of Tables Table 4.1 Indentified parameters for the model (4.1) . . . . . . . . . 56 viii List of Figures Figure 2.1 Viability constraints for the double integrator plotted for k1 = k2 = 0.5. . . . . . . . . . . . . . . . . . . . . . . . . . 12 Figure 2.2 Viability kernel for the double integrator plotted for u0 = 0.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 Figure 2.3 Iteratively constructing an under-approximation of V iab[0,τ ](K). 21 Figure 3.1 A compact, convex set A (white) and the corresponding tight over-approximation (grey) in the set of directions L shown at left. . . . . . . . . . . . . . . . . . . . . . . . . . 30 Figure 3.2 A compact, convex set A (white) and the corresponding tight under-approximation (grey) in the set of directions L shown at left. . . . . . . . . . . . . . . . . . . . . . . . . 31 Figure 3.3 Viability constraints K = {x : ||x||∞ ≤ 0.5} (light blue) and the corresponding under-approximation of the viabil- ity kernel V iab(K) (grey) for the double integrator (3.10) computed using Algorithm 3.2. The approximation is per- formed using a fixed-complexity polytope with 8 vertices. 42 Figure 3.4 Viability constraints K = {x : ||x||2 ≤ 0.5} (light blue) and a corresponding single-ellipsoid under-approximation of the viability kernel V iab(K) (grey) for the double inte- grator (3.10) computed using Algorithm 3.3. The level set approximation (which is a highly accurate approximation of the continuous-time viability kernel, and is computa- tionally feasible for this low-dimensional system) is shown in dark blue. . . . . . . . . . . . . . . . . . . . . . . . . . 43 ix Figure 3.5 Viability constraints K = {x : ||x||2 ≤ 0.5} (light blue) and a corresponding support vector under-approximation of the viability kernel V iab(K) (grey) for the double inte- grator (3.10) computed using Algorithm 3.4. The level set approximation (which is a highly accurate approximation of the continuous-time viability kernel, and is computa- tionally feasible for this low-dimensional system) is shown in dark blue. . . . . . . . . . . . . . . . . . . . . . . . . . 44 Figure 3.6 Viability constraints K = {x : ||x||∞ ≤ 0.5} (light blue) and the corresponding under-approximation of the viabil- ity kernel V iab(K) (grey) for the double integrator (3.10) computed using Algorithm 3.2. The approximation is per- formed using a polytopes with varying numbers of vertices. 45 Figure 3.7 Viability constraints K = {x : ||x||2 ≤ 0.5} (light blue) and a corresponding piecewise-ellipsoidal under-approximation of the viability kernel V iab(K) (grey) for the double inte- grator (3.10) computed using Algorithm 3.3. The level set approximation (which is a highly accurate approximation of the continuous-time viability kernel, and is computa- tionally feasible for this low-dimensional system) is shown in dark blue. . . . . . . . . . . . . . . . . . . . . . . . . . 46 Figure 3.8 Viability constraints K = {x : ||x||2 ≤ 0.5} (light blue) and a corresponding support vector under-approximation of the viability kernel V iab(K) (grey) for the double inte- grator (3.10) computed using Algorithm 3.4. The level set approximation (which is a highly accurate approximation of the continuous-time viability kernel, and is computa- tionally feasible for this low-dimensional system) is shown in dark blue. . . . . . . . . . . . . . . . . . . . . . . . . . 47 Figure 3.9 Comparison of the run time for a chain of integrators of length n. . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 x Figure 3.10 Viability kernel for the heat conduction problem projected onto various coordinate planes. The projection of the state constraint set is shown in red and under- and over- approximations of the viability kernel are shown in blue and grey respectively. This approximation in 120 direc- tions was performed in t = 425.16 seconds. As expected, we see that the set of viable states appears smaller when projected onto coordinates farther from the heated end at ξ1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 Figure 4.1 Transfer diagram for the respiratory system model (4.1). . 55 Figure 4.2 Phase plane for system (4.1) with u(t) = umax. The tra- jectory forming the right-most boundary of the viability kernel is shown in bold. . . . . . . . . . . . . . . . . . . . 57 Figure 4.3 Phase plane for system (4.1) with u(t) = 0. The trajec- tory forming the lower-left-most boundary of the viability kernel is shown in bold. . . . . . . . . . . . . . . . . . . . 58 Figure 4.4 Viability kernel for system (4.1) computed using Algo- rithm 3.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 Figure 4.5 Viability kernel for the 6-dimensional pharmacokinetic model given in Subsection 4.3 computed using Algorithm 3.4. Constraints and viability kernels are shown projected onto the first 3 coordinates of the 6-dimensional state space. These approximations were computed in t = 693s. 61 xi Acknowledgements The work summarized in this thesis has been possible only thanks to the contributions of many different people who, in their different ways, helped, supported and guided me over the course of my Master’s degree program. In particular, I would like to thank my advisor Guy Dumont for his willingness to take a poor young mathematician off the streets and help him to find his way in the strange and frightening world of control engi- neering. Many thanks to Shahab Kaynama for introducing me to viability theory and providing the inspiration for much of the work in this thesis, and for many discussions on everything from first-order logic to optimization to computational geometry. I am grateful to Meeko Oishi for sharing with me her expertise in reachability, and for her guidance in navigating the aca- demic world. Thanks also to Ian Mitchell for providing abundant, insightful commentary on the many iterations of the manuscripts that eventually be- came this thesis. Finally, to all my fellow students in UBC’s Electrical & Computer Engineering in Medicine group: thank you, it has been a jolly time. xii To my parents—for instilling in me a love of learning xiii Chapter 1 Introduction As digital processor technology continues to improve, computational ele- ments are increasingly being embedded in medical devices. Among other benefits, this technology enables computationally intensive control elements (e.g. model predictive controllers) to be incorporated as part of a medi- cal device. Physiologic closed-loop control systems, which include a human patient in the feedback loop, have the potential to improve the safety and reliability of medical treatment. Yet the complexity of such systems makes them difficult to analyse. In order to provide mathematical guarantees that the closed-loop control system is safe, new verification techniques are needed. Computational verification of discrete event systems, or model checking, allows engineers to ensure that a hardware or software design meet a given set of specifications [10]. In the context of control systems, verification typ- ically consists of ensuring that the state of the control system will avoid a subset of the state space designated as undesirable, or remain in a subset designated as desirable. As control systems often involve the interaction between discrete- and continuous event components, there has been sub- stantial interest in techniques for formal verification of control systems in the hybrid systems community [9, 28, 39]. Set-based methods for control system verification, which rely on the com- putation of reachable or invariant sets, have received a substantial amount of research attention in recent years (cf. [6] and references therein). These 1 methods allow engineers to provide guarantees regarding the safety and re- liability of systems because set-based methods allow for the study of all possible behaviours of a system. To demonstrate that a system is safe, we consider the space of all possible system configurations, or states. The con- figuration space is divided into “safe” and “unsafe” states by specifying a set of “unsafe” states that we wish to avoid, or by specifying a set of “safe” states that we wish to remain inside. For system models containing an external input, it is important to distin- guish between two interpretations of the input. Co-operative inputs, such as control inputs, are beneficial to the safety of the system in question. As such, for a given starting configuration we ask whether there exists a co-operative input that will keep the system safe. On the other hand, adversarial inputs, such as noise, disturbances or time-varying uncertainties are detrimental to the system’s safety. In this case we need to ensure that for every adversarial input, the system will remain safe. 1.1 Computational approaches to verification We consider a system modelled by a d-dimensional system of first-order ordinary differential equations ẋ = f(x, u) (1.1) where the d-dimensional state vector x parametrizes the possible configura- tions, or states, of our system and the time-varying external input u takes values in a set U ⊆ Rm. 1.1.1 Simulation The most straight-forward method for verification involves simulating in- dividual trajectories of our system. Any trajectory that leaves the set of safe states (or equivalently enters the set of unsafe states) is an undesirable trajectory. Faced with adversarial inputs, we can demonstrate that a particular 2 starting configuration is unsafe by finding a single input that drives the system into an unsafe configuration. Since it is only possible in practice to simulate a finite number of starting configurations with a finite number of inputs, we cannot use this method guarantee that a system is safe for all initial configurations despite all possible inputs. However, a single trajectory leaving the safe region demonstrates that the system is unsafe. Though it is impossible to guarantee safety via this method, we could provide evidence that the system is likely safe by simulating a large number of trajectories. For co-operative inputs, we can demonstrate that a starting configuration is safe by finding any input that keeps the system in the safe region for all time. Again, it is not possible to explore all possible inputs to the system and starting configuration, but simulation can be used to provide evidence for the existence of a safe input for each starting configuration. 1.1.2 Barrier functions and invariant sets A more rigorous approach to the verification of systems with adversarial inputs involves establishing the existence of a mathematical barrier between the set of possible initial configurations and the set of unsafe states. We can certify that the system is safe by demonstrating that no trajectory can cross the barrier to reach the set of unsafe states. Geometrically, this is interpreted to mean that at each point x on the barrier, the cone of velocities F (x) = {f(x, u) : u ∈ U} lies entirely on the safe side of the barrier. One way of establishing the existence of such a barrier is to find a po- tential function, or barrier certificate [35], V : Rd → R whose zero level set separates the initial states from the unsafe states. The function V should take non-positive values for each initial state and take positive values for each unsafe state. The barrier between these sets is given by the level set {x : V (x) = 0}. For adversarial inputs, for every point x on the barrier and every input u, the time derivative of the barrier function should satisfy V̇ (x) = ∂V ∂x · f(x, u) ≤ 0, (1.2) indicating that all velocity vectors point towards the safe side of the barrier. 3 For co-operative inputs, for each point x on the barrier, there needs to exist some input u such that (1.2) is satisfied, indicating that there is some velocity vector that points towards the safe side of the barrier. The search for a barrier function presents many of the same challenges as the search for a Lyapunov function. For polynomial systems this search can be automated through sum of squares methods [34]. 1.1.3 Reachability-based approaches Rather than simulating individual trajectories of the system (1.1), reachability- based approaches to verification compute a “trajectory of sets” correspond- ing to all possible initial conditions and inputs to the system [6]. This method can be thought of as a breadth first approach to simulation [25]. Its advantage over depth first simulation comes from an increased ability to make guarantees about system safety. In the maximal forward reachability-based approach to verification, we compute the set of all states that can be reached from the initial set of starting configurations. If this set does not intersect the set of unsafe states, we can certify that the system is safe despite any adversarial input that could be applied. We can also use minimal backward reachability techniques. Given a set of unsafe states, we compute the set of starting configurations which reach this unsafe set no matter what co-operative input is applied. If no starting configuration lies in this set, this ensures that for each starting configuration there exists some co-operative input keeping the system within the safe region. 1.2 Viability kernel approximation methods In systems with inputs, the viability kernel is an example of the invariant set approach. Viability theory (cf. [2]; [4]) plays an important role in safety verification for control systems. Consider, for example, problems in aerodynamic flight envelope protection [39] or underwater vehicle operation under constraints [33]. Because the control authority for the system is also 4 constrained, there are some configurations in the safe set for which the state may inevitably exit. Hence it is important to identify the subset of the safe set for which the existence of a control input that keeps the state of the system within the safe region can be guaranteed. This subset, known as the viability kernel (or as the controlled invariant set), takes into account the system’s dynamics and bounded control authority. For a constraint set K, the viability kernel V iab(K) is the subset of K that consists of all states for which a control input exists that keeps the state of the system within K for the duration of a known (possibly infinite) time horizon. The viability kernel has traditionally been computed using Eulerian methods such as the Viability Kernel Algorithm [38] and level set approaches [32]. However, Eulerian methods involve numerical implementations that re- quire gridding the state space and thus their time and memory complexity grow exponentially with the state dimension. In practice, this approach is infeasible for systems with more than 3 or 4 states. There has been recent work to approximate the viability kernel for high-dimensional systems. We survey various viability kernel approximation methods in this section. 1.2.1 Viability kernel algorithm In [38], Saint-Pierre presents an algorithm for approximating the viability kernel recursively as a sequence of finite-horizon viability kernels Kn of in- creasing horizon length n. He also presents a thorough discussion of the algorithm’s convergence properties. The details of the viability kernel algo- rithm are given in Section 2.2. 1.2.2 Supervised classification In [12], Deffuant, Chapel and Martin describe a method of computing the viability kernel in systems with a high-dimensional input space. The method uses a supervised classification procedure (e.g. a support vector machine) associating a finite set of training tuples in K × {±1} with a classification function ` : K → {±1} to enable the approximate computation of the convex hull L(Khn+1) of a set K h n+1 = {x ∈ Khn : d(G(x), L(Khn)) < } of grid points 5 whose convex hull approximates the Saint-Pierre set Kn+1. This method enables the analysis of systems with a high-dimensional input space but since it still relies on gridding the state space, its scalability in terms of state dimension is still limited. In [8] this method is applied to the analysis of an ecological system with a 51-dimensional input but only 6 state dimensions. 1.2.3 Approximate dynamic programming The viability kernel is computed in [11] as the zero level set of the viscosity solution V ∗ of the Hamilton-Jacobi-Bellman equation V ∗(x) log(γ) + min u∈U [∇V ∗(x) · f(x, u)] = 0. (1.3) with boundary condition V ∗(x) = 1 for x 6∈ K, and where γ ∈ (0, 1). This equation is discretized with time step ∆t and the approximate solution V ∆t of the resulting Bellman equation V ∆t(x) = γ∆t min u∈U V ∆t(x+ f(x, u)∆t) is computed via approximate dynamic programming. Given > 0, an ap- proximation to the viability kernel is given by the sublevel set K = {x ∈ K : V ∆t(x) ≤ }. In [11], this algorithm is applied to a four-dimensional anaerobic diges- tion model. Using a grid of 2×105 points sampled uniformly in the viability kernel, the viability kernel takes 162.52 sec. to compute. 1.2.4 Simulated annealing Bonneuil [7] formulates the viability constraint set K as the zero sublevel set of a function h: K = {x : h(x) ≤ 0}. If S(x) denotes the set of trajectories of the system originating at x, the viability kernel over a time horizon τ is given by V iab(K) = { x : inf x(·)∈S(x) sup t∈[0,τ ] h(x(t)) ≤ 0 } . 6 Thus this problem can be solved by minimizing c(x) = infx(·)∈S(x) supt∈[0,τ ] h(x(t)) by searching S(x) at discrete time intervals with a procedure such as simu- lated annealing. The simulated annealing method has been demonstrated for a chain of integrators in 10 dimensions discretized with 36 steps. The method takes on average 22 minutes to compute each initial point on viability kernel’s boundary. 1.3 Contributions of this thesis This thesis presents a connection between the viability kernel and reachable sets that allows the large class of Lagrangian methods developed for the approximation of maximal reachable sets (e.g. [9, 20, 25]) to be applied for the approximation of viability kernels. In contrast to Eulerian methods, which compute the flow of a vector field based on a coordinate system that is fixed in space, Lagrangian methods use representations that follow the vector field’s flow. Since Lagrangian methods do not depend on gridding the state space, they can be used to analyse high-dimensional systems. We provide four examples of such algorithms. The exact polytope method given in Algorithm 3.1 performs well in terms of accuracy but does not scale well as the state dimension grows, becoming infeasible in greater than four dimensions. In comparison, the time complexity of the approximate polytopic method described in Algorithm 3.2 and the ellipsoidal method described in Algorithm 3.3 increase more slowly with the state dimension, but their accuracy is limited due to the under-approximations made. The support vector method given in Algorithm 3.4 strikes a balance between scalability and accuracy. It allows the user to choose a desired accuracy in terms of the number of points on the boundary of the viability kernel that they wish to evaluate. Near the end of this thesis, we discuss an application of viability meth- ods to problems arising in the closed-loop control of physiological systems. We argue that the viability kernel can be used to detect safety violations in advance, allowing a fallback mode of operation to be initiated early. We 7 provide two examples of such systems: the control of the carbon dioxide partial pressure in the human respiratory system under mechanical ventila- tion and the control of the concentration of the anaesthetic drug Propofol in anaesthetized children. 8 Chapter 2 Establishing connections between viability and reachability There is a close relationship between viability theory [4] and constrained reachability [22]. Both frameworks study the evolution of dynamic systems under input and/or state constraints. The relationship between the two the- ories is often discussed in the context of optimal control theory by formu- lating both reachability and viability problems in terms of Hamilton-Jacobi equations, for example [27]. The Hamilton-Jacobi approach has proven extremely successful in the analysis of low-dimensional systems. Level set methods can be used to ap- proximate the viscosity solution of the Hamilton-Jacobi PDE corresponding to a given viability or reachability problem [39] and tools are available for computing viable and reachable sets numerically [31], however they scale poorly with state dimension. The recent emergence of accurate and scalable methods and tools for ap- proximating reachable sets in high-dimensional systems [13, 20] has inspired us to attempt to find analogous methods for the approximation of viability kernels. In this chapter, we expose a new connection between viability the- ory and reachability theory. This characterization leads to new algorithms 9 for the computation of viability kernels for discrete- and continuous-time systems. The results in this section appeared in preliminary form in a con- ference paper [19]. 2.1 Preliminaries We are concerned with analysing systems of the form{ L(x(t)) = f(x(t), u(t)) u(t) ∈ U (2.1) where the time t ranges throughout a time scale T. The time scale T can be either continuous (T = [0, τ ] ⊆ R+) or discrete (T = [0, τ ] ∩ Z+). If 0 < τ < ∞ this problem is said to have a finite time horizon where as if τ = ∞, it is said to have an infinite time horizon. L is the differential operator corresponding to the given time scale (differentiation in the case of a continuous-time system and first differencing in the case of a discrete-time system). The system’s state x ranges over the finite-dimensional vector space Rn and the system’s input is constrained to a nonempty, compact, convex subset U ⊆ Rm.1 When Equation (2.1) evolves under continuous time, we assume that the function f : Rn × U → Rn is sufficiently smooth to guarantee the existence and uniqueness of solutions to the corresponding initial value problem. It is often convenient to formulate problem (2.1) in an equivalent way using the language of set-valued maps and differential inclusions [3]. We define the set-valued map F : Rn Rn as F (x) = {f(x, u)|u ∈ U}. This allow us to write (2.1) as a differential inclusion L(x(t)) ∈ F (x(t)). (2.2) 1This is not the most general context in which viability theory can be developed. Aubin [2] allows the constraint set U to be dependent on the state x. 10 2.1.1 Elements of viability theory Viability theory is concerned with ensuring that a system’s state x remains within a set of viability constraints K ⊆ Rn. Essentially, any trajectory of system (2.1) that leaves the set K at some point in time is considered to be no longer viable. We call a set S viable under system (2.1) if for every initial state x0 ∈ S there exists some measurable input u0 : T → U such that the solution x(·) to the initial value problem{ L(x(t)) = f(x(t), u0(t)) x(0) = x0 (2.3) satisfies x(t) ∈ S for all t ∈ T. The viability kernel of a set of viability constraints K is the largest viable set contained in K. Equivalently, the viability kernel is defined as follows: Definition 1. V iabT(K) = {x0 ∈ K | ∃u0 : T→ U ∀t ∈ T x(t) ∈ K}. Essentially, the viability kernel provides a refinement of the set of via- bility constraints that takes the systems dynamics (2.1) into account. This interpretation is illustrated in Example 1. Example 1. Consider a point mass travelling in one dimension. Let the state x1 ∈ R denote its position and x2 ∈ R its velocity. We consider a state to be viable (or “safe”) if |x1| ≤ k1 and |x2| ≤ k2. The corresponding set of viability constraints K is shown in Figure 2.1. We assume that the mass’ dynamics are that of a double integrator and that only a finite force can be applied to it and hence its acceleration is constrained to have magnitude less than u0. These dynamics are described 11 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 x1 x 2 Figure 2.1: Viability constraints for the double integrator plotted for k1 = k2 = 0.5. by Equation (2.4): [ ẋ1(t) ẋ2(t) ] = [ 0 1 0 0 ][ x1(t) x2(t) ] + [ 0 1 ] u(t) u(t) ∈ U = [−u0, u0]. (2.4) The dynamics play an important role in determining which states will remain safe as the system evolves over time. For example, if we begin a trajectory with the mass moving in the positive direction at velocity v, it will travel a distance of at least d(v) = 1 2u0 v2 before stopping regardless of the input applied. Therefore, to be able to guar- antee the system’s safety, we require that the mass’ starting position x1(0) be at least a distance d(v) from the border of the safe set at x1 = k1. That is, we require that k1 − x1(0) ≥ 1 2u0 v2. 12 Applying the same argument for a negative initial velocity v, we get the additional requirement k1 − x1(0) ≥ 1 2u0 v2. The set of states satisfying these additional safety constraints make up the viability kernel (shown for this system in Figure 2.2). The viability kernel consists of the set of states that can be made to remain safe throughout the time scale T by choosing a an appropriate input u : T→ U . −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 x1 x 2 Figure 2.2: Viability kernel for the double integrator plotted for u0 = 0.3. 2.1.2 Reachability under input constraints Constrained reachability analysis is a common technique for formal safety verification (for example 30). It provides a method of simulating all possible trajectories of a dynamic system under all admissible inputs. Essentially, it is concerned with determining if any trajectories of the system (2.1) that begin in a set of initial conditions I can reach a set of terminal states T . There are two ways to approach the problem of reachability analysis. We can begin by considering the set of initial states and follow this set forward in time under the flow of (2.1) to compute what is known as the forward 13 reachable set. The other approach considers the set of terminal states T and follows the flow of (2.1) backward in time to compute the backward reachable set. For our purposes, it is appropriate to use the backward approach. We define the set backward reachable from T over a time domain T with finite horizon τ as follows: Definition 2. ReachT(T ) = {x0 ∈ Rn | ∃u0 : T→ U x(τ) ∈ T} where x(·) denotes the solution to the initial value problem (2.3). 2.2 The Viability Kernel Algorithm We are interested in computing viability kernels numerically. To do this, we require a set-valued numerical algorithm that can be made to scale well with the dimension of the state space. The classical Viability Kernel Al- gorithm due to Patrick Saint-Pierre [38] provides a stepping stone towards such an algorithm. It can be used to compute the finite-time viability kernel for discrete-time systems as well as an approximation for continuous-time systems by using a discretization of the continuous dynamics. As it will form the basis for our algorithm, we review Saint-Pierre’s algorithm in this section and establish some of its properties. We consider a discrete-time differential inclusion with time step ρ > 0 x(t+ 1)− x(t) ρ ∈ F (x(t)). (2.5) By defining the new set-valued map Gρ : Rn Rn Gρ(x) = {x+ ρF | F ∈ F (x)} we can write Equation (2.5) as x(t+ 1) ∈ Gρ(x(t)). (2.6) 14 Saint-Pierre computes the viability kernel of a compact set K under (2.6) as follows. We define a sequence of sets recursively by{ K0 = K Kn+1 = {x ∈ Kn | Kn ∩Gρ(x) 6= ∅}. (2.7) This defines a decreasing sequence of sets K0 ⊇ K1 ⊇ K2 ⊃ · · · . Proposition 1 characterizes these sets as finite-horizon viability kernels. Proposition 1 (Characterization of Saint-Pierre sets as finite-horizon via- bility kernels). The τ -th set Kτ defined in Equations (2.10) represents the viability kernel with time horizon τ . That is, Kτ = V iabT(K) where T = Z ∩ [0, τ ]. Proof. We prove Proposition 1 by induction on τ . To establish the base, note that by definition V iab{0}(K) = K = K0. Now suppose that for some τ ∈ Z+ we have V iab[0,τ ]∩Z(K) = Kτ . We wish to show that V iab[0,τ+1]∩Z(K) = Kτ+1. Let x0 ∈ V iab[0,τ+1]∩Z(K). Then there exists an input u0 : [0, τ+1]∩Z→ U such that the solution x(·) to the initial value problem (2.3) satisfies x(t) ∈ K for all t ∈ [0, τ + 1] ∩ Z. Define x1 = x(1) = x0 + ρf(x0, u0(0)). By definition, x1 ∈ Gρ(x0). Further, since we know that x(·) solves the IVP (2.3), we have x1 ∈ V iab[1,τ+1]∩Z(K) = Kτ . Thus Kτ ∩Gρ(x0) 6= ∅ so x0 ∈ Kτ+1. Conversely, let x0 ∈ Kτ+1. Then Kτ ∩Gρ(x0) 6= ∅, or equivalently, there is some point x1 with x1 ∈ Kτ and x1 = x0 + ρf(x0, u∗) for some u∗ ∈ U . Since x1 ∈ Kτ = V iab[0,τ ]∩Z(K), there exists an input u1 : [0, τ ] ∩ Z → U such that x(t) ∈ K for t ∈ [0, τ ] ∩ Z. We define u0(t) = { u∗ if t = 0 u1(t− 1) if t ∈ [1, τ + 1] ∩ Z Using this input, the solution x(·) to the IVP (2.3) satisfies x(t) ∈ K for all 15 t ∈ [0, τ + 1] ∩ Z. So x0 ∈ V iab[0,τ+1]∩Z(K). We can also make use of the sets Kn to approximate the infinite-horizon viability kernel. We define K∞ = ⋂∞ n=0Kn. The following theorem is due to Saint-Pierre: Theorem 1 (Saint-Pierre, [38]). Let K be a compact set and let V iabGρ(K) be the infinite-horizon viability kernel of K under (2.6). Then K∞ = V iabGρ(K). Thus the sequence of sets Kn converges from above to the infinite-horizon viability kernel. Thus we can compute an approximation to the infinite- horizon viability kernel by computing Kn for n sufficiently large. Further- more, the discrete-time viability kernel that we compute here using the explicit discretization (2.5) under-approximates the viability kernel of the original continuous-time system ẋ(t) ∈ F (x(t)). (2.8) As the discretization time step ρ → 0, the approximations converge to a viable subset of the continuous-time viability kernel in a sense made precise by Theorem 2. Theorem 2 (Saint-Pierre, [38]). The upper limit lim sup ρ→0 V iabGρ(K) is a viable set under F . That is, lim sup ρ→0 V iabGρ(K) ⊆ V iabF (K). Theorem 2 guarantees that we can compute an under-approximation of the viability kernel for a continuous-time system by using a discretized algorithm with a sufficiently small time step. 16 2.3 Computing viability kernels using reachability techniques In this section, we present a method of expressing finite horizon viability kernels in terms of reachable sets. This provides a modified version of Saint- Pierre’s Viability Kernel Algorithm that can be implemented using efficient and scalable techniques developed within the context of reachability analysis. We present the algorithm first for discrete-time systems then for continuous- time systems. In the discrete case, the viability kernel is computed exactly, while in the continuous-time case, we compute an under-approximation of the true viability kernel. Throughout this thesis, whenever we cannot compute a set exactly we ensure that our approximation is an under-approximation. This guarantees that the set we compute to approximate V iab(K) is indeed viable under K. Discrete-time systems We consider the case when the system (2.1) evolves in discrete time. The system’s dynamics can be described by the constrained difference equation{ x(t+ 1) = f(x(t), u(t)) u(t) ∈ U . (2.9) For discrete-time systems, the viability kernel can be computed using Saint- Pierre’s Viability Kernel Algorithm via the following recursive formula that gives the finite-horizon viability kernel Kn+1 := V iab[0,n+1]∩Z+(K) in terms of the finite-horizon viability kernel Kn at the previous step [38]:{ K0 = K Kn+1 = {x ∈ Kn | Kn ∩ F (x) 6= ∅} (2.10) where F (x) = {f(x, u) | u ∈ U}. We can reformulate this recursive definition of the finite horizon viability kernels Kn in terms of reachable sets with time horizon τ = 1: Theorem 3. The sequence of finite horizon viability kernels Kn can be 17 computed recursively in terms of reach sets as{ K0 = K Kn+1 = K0 ∩Reach1(Kn). (2.11) Proof. Let the constrained difference equation (2.9) be expressed as the difference inclusion x(t+ 1) ∈ F (x(t)). Then using the definition of Kn+1, x ∈ Kn+1 ⇐⇒ x ∈ Kn ∧ (F (x) ∩Kn 6= ∅) ⇐⇒ x ∈ Kn ∧ ∃y(y ∈ F (x) ∧ y ∈ Kn) ⇐⇒ x ∈ Kn ∧ ∃y(∃u ∈ U y = f(x, u) ∧ y ∈ Kn) ⇐⇒ x ∈ Kn ∧ ∃u ∈ U f(x, u) ∈ Kn ⇐⇒ x ∈ Kn ∧ x ∈ Reach1(Kn) ⇐⇒ x ∈ Kn ∩Reach1(Kn). Thus Kn+1 = Kn ∩ Reach1(Kn). We prove that Kn ∩ Reach1(Kn) = K0 ∩ Reach1(Kn) by induction. The base K0 ∩Reach1(K0) = K0 ∩Reach1(K0) is clear. It follows from Kn+1 ⊆ Kn that K0 ∩ Reach1(Kn+1) ⊆ K0 ∩ Reach1(Kn) = Kn∩Reach1(Kn) = Kn+1. Hence we have the first inclusion K0∩Reach1(Kn+1) ⊆ Kn+1∩Reach1(Kn+1). The opposite inclusion follows from the fact that Kn+1 ⊆ K0. This recursive formula leads to Algorithm 2.1 for computing the finite horizon viability kernel over the discrete time interval T = {t ∈ Z+ | t ≤ N}. 18 Algorithm 2.1 Exact computation of the viability kernel (discrete-time) K0 ← K n← 0 while n ≤ N do if Kn = ∅ then . If true, V iabT(K) = ∅ KN ← ∅ break end if if Kn = Kn−1 then . If true, V iabT(K) = Kn KN ← Kn break end if L← Reach1(Kn) Kn+1 ← K0 ∩ L n← n+ 1 end while return (KN ) . KN = V iabT(K) Continuous-time systems We now consider the case when the system (2.1) evolves in continuous time. In this case, the system’s dynamics are described by the constrained differ- ential equation { ẋ(t) = f(x(t), u(t)) u(t) ∈ U . (2.12) Before we can present our algorithm, we require a few definitions. We say that a vector field f : Rd × U → Rd is bounded on K ⊆ Rd if there exists a norm || · || : Rd → R+ and a real number M > 0 such that for all x ∈ K and u ∈ U we have ||f(x, u)|| ≤ M . We also define the || · ||-distance of a point x ∈ Rd from a nonempty set S ⊂ Rd as dist(x, S) = inf s∈S ||x− s||. (2.13) 19 Computing an under-approximation of the viability kernel Let K be the set of viability constraints for the constrained differential equation (2.12). We assume that the vector field f is bounded by M on K in the norm || · ||. Given a discretization time interval ρ, we begin by defining an under-approximation of the viability constraint set (Figure 2.3a): Kρ := {x ∈ K | dist(x,Kc) ≥ ρM}. (2.14) We under-approximate K by a distance ρM because we are only considering the system’s state at discrete times tn = nρ. At a time t in the interval [tn, tn+1], a solution x(·) of (2.12) can travel a distance of at most ||x(tn)− x(t)|| ≤ ∫ t tn ||ẋ(τ)||dτ ≤M(t− tn) ≤ ρM from its initial location x(tn). Thus the under-approximation (2.14) ensures that the state does not leave K at any time in the interval [tn, tn+1]. We proceed by defining a sequence of sets Kn(ρ) analogously to the discrete-time case. The under-approximation Kρ of the viability constraints defines the base of our recursion. We then define subsequent sets Kn(ρ) recursively: { K0(ρ) = Kρ Kn+1(ρ) = K0(ρ) ∩Reach[0,ρ](Kn(ρ)). (2.15) At each time step, we calculate the set of states from which Kn(ρ) is reach- able, then intersect this set with the set of safe states. This process is illustrated in Figure 2.3. Each computed set Kn(ρ) is an approximation of the finite horizon viability kernel V iab[0,τ ](K) for τ = nρ. Note that the resulting set depends on our choice of the time step ρ. We claim that for any ρ > 0, Kn(ρ) under-approximates V iab[0,nρ](K). Theorem 4. Suppose that the vector field f : Rd×U → Rd is bounded on a set K ⊆ Rd. Then for any time step ρ the sets Kn defined by the recurrence 20 (a) We define the initial under- approximation of the safe set K0(ρ) = Kρ. (b) We calculate the set of backward reach- able states from K0(ρ). (c) We intersect the backward reachable set with the initial set to get K1(ρ) (d) Next, we calcu- late the set of back- ward reachable states from K1(ρ). (e) Again, we inter- sect the backward reach- able set with the ini- tial set to get a new set K2(ρ) (f) By repeat- ing this process, we eventually reach an under-approximation Kn(ρ) of the viability kernel. Figure 2.3: Iteratively constructing an under-approximation of V iab[0,τ ](K). 21 relations (2.15) satisfy Kn(ρ) ⊆ V iab[0,nρ](K). (2.16) Proof. Since f is bounded on K, there exists a norm || · || and a real number M > 0 with ||f(x, u)|| ≤ M for all x ∈ K. Now, take a point x0 ∈ Kn(ρ). By the construction of Kn(ρ), this means that for each k = 1, . . . , n there is some point xk ∈ Kk(ρ) and an input uk : [0, ρ] → U such that xk can be reached from xk−1 at time ρ using input uk. Thus, taking the concatenation of the inputs uk, we get an input u : [0, nρ]→ U such that the solution x : [0, nρ]→ Rd to the initial value problem ẋ = f(x, u), x(0) = x0, satisfies x(kρ) = xk ∈ Kk(ρ) ⊆ {x ∈ K | dist(x,Kc) ≥Mρ}. We claim that this guarantees that x(t) ∈ K for all t ∈ [0, nρ]. Indeed, any t ∈ [0, nρ) lies in some interval [tk, tk+1) = [kρ, (k + 1)ρ). Since f is bounded by M , we have ||x(tk)− x(t)|| ≤ ∫ t tk ||ẋ(τ)||dτ ≤M(t− tk) ≤ ρM Further, x(tk) ∈ Kk(ρ) implies that dist(x(tk),Kc) ≥ ρM . Combining these, we see that dist(x(t),Kc) ≥ dist(x(tk),Kc)− ||x(t)− x(tk)|| > ρM − ρM = 0 and hence x(t) ∈ K. Thus, x0 ∈ V iab[0,nρ](K). Accuracy of the Approximation By choosing a sufficiently small time step, the approximation can be made arbitrarily accurate in the following sense: fix a time horizon τ and partition the interval [0, τ ] into N subintervals [tn−1, tn] of width ρN = τ/N . The union of the approximating sets KN (ρN ) taken over all N ∈ N is bounded between the viability kernels of K and its interior K̊. Theorem 5. Suppose that the vector field f : Rd × U → Rd is bounded on 22 a set K ⊆ Rd. Then we have V iab[0,τ ](K̊) ⊆ ⋃ N∈N KN (ρN ) ⊆ V iab[0,τ ](K). (2.17) Proof. The second inclusion ⋃ N∈NKN (ρN ) ⊆ V iab[0,τ ](K) follows di- rectly from Theorem 4. To prove the first inclusion, take x0 ∈ V iab[0,τ ](K̊). Then there exists an input u : [0, τ ] → U such that the solution x(·) to the initial value problem ẋ = f(x, u), x(0) = x0, satisfies x(t) ∈ K̊ for all t ∈ [0, τ ]. Since K̊ is open, for any x ∈ K̊ we have dist(x,Kc) > 0. Further, x : [0, τ ] → Rd is continuous so the function t 7→ dist(x(t),Kc) is continu- ous on the compact set [0, τ ]. Thus, we can define d > 0 to be its minimum value. Now take N large enough such that ρN < d/M . We need to show that x0 ∈ KN (ρN ). First note that N is chosen such that dist(x(t),Kc) > ρNM for all t ∈ [0, τ ]. Hence x(tN−n) ∈ Kn(ρN ) for all n = 0. To show that x(tn−1) ∈ Reach[0,ρN ](Kn(ρN )) for all n = 1, . . . , N , consider the sequence of inputs un : [0, ρN ]→ U defined as un(t) = u(tn−1 + t). (2.18) It is easy to verify that for all n, we can reach x(tn) from x(tn−1) at time tn using input un. Thus, in particular, we have x0 = x(0) ∈ Reach[0,ρN ](KN−1(ρN )). So x0 ∈ KN (ρN ). Hence V iab[0,τ ](K̊) ⊆ ⋃ N∈NKN (ρN ). Corollary 1. When K is open,⋃ N∈N KN (ρN ) = V iab[0,τ ](K). (2.19) Algorithm 2.2 computes an approximation of the viability kernel in the continuous-time case using the recursive formula (2.15). Theorem 4 guar- antees that the computed set always under-approximates the true viability kernel while Theorem 5 guarantees that the approximation is asymptotically tight as the time step ρ→ 0. 23 Algorithm 2.2 Under-approximation of the viability kernel (continuous- time) Choose ρ > 0 . Determines approximation accuracy N ← τ/ρ . Number of time steps K0 ← Kρ . Initial under-approximation of K0 n← 0 while n ≤ N do if Kn = ∅ then . If true, V iab[0,τ ](K) = ∅ KN ← ∅ break end if if Kn = Kn−1 then . If true, V iab[0,τ ](K) = Kn KN ← Kn break end if L← Reach[0,ρN ](Kn) Kn+1 ← K0 ∩ L n← n+ 1 end while return (KN ) . KN ⊆ V iab[0,τ ](K) 2.4 Summary In this chapter, we have presented a method of computing the viability kernel of a discrete-time control system using backward reachable sets. For the continuous-time case, we present an approximate method that ensures that the computed set under-approximates the continuous-time viability kernel. These methods are described in Algorithms 2.1 and 2.2 respectively. Since there are a number of established Lagrangian algorithms for com- puting reachable sets, these results will allow us to develop Lagrangian al- gorithms for computing viability kernels. 24 Chapter 3 Lagrangian algorithms for approximating viability kernels in linear systems In Chapter 2 we demonstrated that the viability kernel of an input-constrained dynamic system can be computed in terms of reachable sets. In this chapter we use this result to develop efficient algorithms for computing or approxi- mating the viability kernel in high-dimensional linear systems. As the viability kernel is often used for safety verification, it is desirable that any approximations made are conservative so that the computed safe set under-approximates the true viability kernel. This way, we can guarantee that any point within the computed kernel is truly a safe initial state. There are numerous Lagrangian algorithms for approximating reachable sets in linear systems. These algorithms rely on particular geometric repre- sentations such as polytopes [9], ellipsoids [21], zonotopes [14], or support functions [26]. Each representation has advantages and disadvantages in terms of representation size, approximation fidelity and ease of performing geometric operations. We begin by surveying and comparing various set representations. We then develop a number of practical implementations of Algorithms 2.1 and 2.2 and compare their performance on a pair of benchmark examples. 25 3.1 Set representations for linear reachability We now restrict our attention to linear dynamics1{ L(x(t)) = Ax(t)− v(t) v(t) ∈ V. (3.1) In the discrete-time case, the backward reachable set over a single time step is computed as Reach(K) = A−1(K ⊕ V). (3.2) Here A−1(·) denotes the preimage of a set under the map A : Rd → Rd. Throughout the remainder of this section, we will assume that A is non- singular, and thus the preimage of A can be calculated simply by applying the linear transformation A−1 to the set K ⊕ V. This is a fair assumption because we are mainly concerned with discrete-time systems that arise from the discretization of continuous time systems. Such systems have a dynamics matrix of the form A = exp(ρAc) which is always invertible. As the operations performed on sets include Minkowski summation, lin- ear transformation and intersection, the ideal set representation would be a class of objects closed under these three operations. We also hope for a rep- resentation under which all three operations can be performed accurately, efficiently and using a constant amount of memory. 3.1.1 Convex polytopes A convex polytope P ⊂ Rd (hereafter simply polytope) is a bounded geo- metric object with flat sides. A polytope can be defined as the convex hull of a finite set of points v1, · · · , vk. This is known as the vertex representation and in this case, P is called a V-polytope. A polytope can also be defined as an intersection of half-spaces P = f⋂ i=1 {x ∈ Rd | ai · x ≤ bi} = {x ∈ Rd | Ax ≤ b} 1A general linear system L(x(t)) = Ax(t) + Bu(t) u(t) ∈ U ⊆ Rm can be written in this form by setting V = −BU = {−Bu|u ∈ U} ⊆ Rd. 26 This is known as the facet representation and in this case, P is known as an H-polytope. The problem of switching between the H-representation and the V- representation is known as the vertex/facet enumeration problem [5]. It is a well-studied problem and many algorithms have been proposed to solve it. These algorithms tend to be slow however, taking O(dkf) time where k is the number of vertices and f the number of facets in the polytope. Polytopes are a good choice of set representation for our purposes be- cause the class of polytopes is closed under all three operations that we wish to perform. Let P1 = {x ∈ Rd | A1x ≤ b1} and P2 = {x ∈ Rd | A2x ≤ b2} be two polytopes described by their H-representation. We can compute the image of P1 under the linear transformation T as TP1 = {x ∈ Rd | A1T−1x ≤ b1}. The intersection of P1 and P2 can also be easily computed simply by com- bining the constraints as follows P1 ∩ P2 = {x ∈ Rd | A1x ≤ b1 and A2x ≤ b2} = { x ∈ Rd : [ A1 A2 ] x ≤ [ b1 b2 ]} . The Minkowski sum however is difficult to compute using theH-representation and is typically done by first converting to the V-representation. Given two polytopes P1 = Conv{v1, · · · , vk} and P2 = Conv{u1, · · · , ul} their Minkowski sum is computed as P1 ⊕ P2 = Conv{vi + uj | i = 1, · · · , k j = 1, · · · , l}. The linear transformation of P1 under T is also easily computed as TP1 = Conv{Tv1, · · · , T vk}. The disadvantage of the V-representation is that it is difficult to perform 27 intersections. The Multi-Parametric Toolbox [24] provides a comprehensive set of tools for computations using polytopes. 3.1.2 Ellipsoids An ellipsoid E ⊂ Rd is a smooth geometric object contained by a bounded quadratic surface. Any ellipsoid E ⊂ Rd can be expressed as the image of the Euclidean unit ball under an affine transformation T : Rd → Rd. E is called nondegenerate if the transformation T is nondegenerate (i.e. invertible). Ellipsoids can also be defined uniquely by a symmetric, positive definite shape matrix Q and a centre q ∈ Rd as E(q,Q) = {x ∈ Rd | (x− q) ·Q−1(x− q) ≤ 1}. The class of ellipsoids is closed under nondegenerate linear transformations but it is not closed under either Minkowski summation or intersection. How- ever, the Ellipsoidal Toolbox [23] provides a set of routines to efficiently compute tight under- and over-approximations of Minkowski sums and in- tersections of ellipsoids. 3.1.3 Support functions An arbitrary compact, convex set A ⊂ Rd can be represented in terms of its support function. The support function ρA : Rd → R of A is a convex function defined as ρA(`) = max x∈A x · `. (3.3) The support function ρA is a complete representation of A in the sense that A can be reconstructed from ρA as the intersection of all its supporting half spaces A = ⋂ `∈Rd {x ∈ Rd | x · l ≤ ρA(`)}. (3.4) Support functions are convenient for our purposes because all three oper- ations that we wish to perform can be performed directly on the support 28 functions. This fact is given in Theorem 6. Theorem 6. Let A ⊂ Rd and B ⊂ Rd be compact, convex sets and let A : Rd → Rd be a linear transformation represented by a matrix A. We have the following properties: • ρAB(`) = ρB(AT `) • ρA⊕B(`) = ρA(`) + ρB(`) • ρA∩B(`) = infw∈Rd{ρA(`− w) + ρB(w)}. Proof. The first two properties follow directly from the definition. The proof of the third is given in [37]. Thus given support functions for A and B, it is simple to compute the support functions of A⊕ B, AB and A ∩ B. A compact, convex set A ⊂ Rd can be over-approximated by an H- polytope with arbitrary accuracy by sampling its support function. Con- sider a finite set of vectors L ⊂ Rd. Following [25], we can define an over- approximation of A by restricting the intersection in equation (3.4) to the set L, giving us A↑ = ⋂ `∈L {x ∈ Rd | x · ` ≤ ρA(`)}. This over-approximation is tight (i.e. the approximation touches the bound- ary of A) in the directions of L. An example of this approximation is shown in Figure 3.1. 3.1.4 Support vectors When computing viability kernels, it is usually desirable to under-approximate rather than over-approximate a given set. Again following [25], for a com- pact, convex set A ⊂ Rd we can construct an under-approximation of A using support vectors. Given a direction vector ` ∈ Rd the set of support vectors of A is defined as vA(`) = argmax x∈A x · `. (3.5) 29 Figure 3.1: A compact, convex set A (white) and the corresponding tight over-approximation (grey) in the set of directions L shown at left. The support vectors and support function of a convex set are related via the subgradient operation ∂ (see Appendix) as vA(`) = ∂ρA(`). In particular, when ρA is differentiable at `, the set of support vectors in the direction ` is the singleton set vA(`) = {∇ρA(`)}. As was the case for support functions, all three operations that we wish to perform can be performed directly on the support vectors. This result is given in Theorem 7. Theorem 7. Let A ⊂ Rd and B ⊂ Rd be compact, convex sets and let A : Rd → Rd be a linear transformation represented by a matrix A. We have the following properties: • vAA(`) = AvA(AT `) • vA⊕B(`) = vA(`)⊕ vB(`) • vA∩B(`) = vA(`− w̄) ∩ vB(w̄) where w̄ = arg infw∈Rd{ρA(`− w) + ρB(w)}. Proof. Again, the first two properties follow directly from the definition. The proof of the third requires some background in convex analysis and is given in the Appendix. 30 Figure 3.2: A compact, convex set A (white) and the correspond- ing tight under-approximation (grey) in the set of directions L shown at left. As was the case for support functions, the set A can be reconstructed from its support vectors: A = Conv ⋃ `∈Rd vA(`) . Now, given a subset L of directions, we can define an under-approximation of A that is tight in the directions L by selecting a support vector u` ∈ vA(`) in each direction ` ∈ L. An under-approximation of A is then given by A↓ = Conv({u` | ` ∈ L}). This approximation is illustrated in Figure 3.2. 3.2 Algorithms We now present three algorithms that compute an approximation of the viability kernel of the discrete-time system{ x(t+ 1) = Ax(t)− v(t) v(t) ∈ V. (3.6) 31 We compare these three algorithms in Section 3.3. For continuous-time systems, an under-approximation of the viability kernel can be computed using these algorithms, provided that we have an appropriate discretization of the system’s dynamics. 3.2.1 Exact polytopic method The class of polytopes is closed under linear transformation, Minkowski summation and intersection. Therefore when the input constraint set V and the viability constraints K are both polytopes, we can compute the viability kernel of the system (3.6) exactly. This method is similar to other methods of computing viability kernels and controlled invariant sets as described in [6] and implemented in the Multi-Parametric Toolbox [24]. Algorithm 3.1 computes the viability kernel using Algorithm 2.1 with Reach(Kn) computed using (3.2). Algorithm 3.1 Exact polytopic method K0 ← K n← 0 while n ≤ N do if Kn = ∅ then KN ← ∅ break end if if Kn = Kn−1 then KN ← Kn break end if L← A−1(Kn ⊕ V) Kn+1 ← K0 ∩ L n← n+ 1 end while return (KN ) . KN = V iabT(K) 32 Since no approximations are made, the accuracy of this algorithm is per- fect. However, the amount of information required to represent the polytope Kn increases exponentially with successive Minkowski sums as the number of vertices of Kn ⊕ U is (in the worst case) |V (Kn)| · |V (U)|. 3.2.2 Approximate polytopic method Algorithm 3.1 can be used to compute the viability kernel of a discrete-time system exactly using polytopes. However, its run time increases exponen- tially with the time horizon because the number of vertices in the polytopes Kn increases exponentially. We can circumvent this problem by introducing approximations of the intersection and Minkowski sum operations that keep the number of vertices constant. Under-approximating the Minkowski sum of two sets Let K be a convex polytope represented as the convex hull of its vertices K = Conv({ki|i ∈ I}) and let V be a compact, convex set. We wish to find an under-approximation of K ⊕ V that is no more complex that K. The following algorithm calculates an under-approximation Approx⊕(K,V) of K ⊕ V that is expressed as the convex hull of |I| points. 1. Find the centroid (or centre of mass) k∗ of K. 2. For each vertex ki of K: (a) find a support vector vi ∈ vV(ki−k∗) that is tight in the direction ki − k∗. (b) let k̃i = ki + vi. 3. Let Approx⊕(K,V) = Conv({k̃i|i ∈ I}). Under-approximating the intersection of two close polytopes Let K = Conv({ki|i ∈ I}) and K ′ = Conv({k′i|i ∈ I}) be two convex poly- topes. 33 Definition 3. We say that two convex polytopes K = Conv({ki|i ∈ I}) and K ′ = Conv({k′i|i ∈ I}) are close if max i∈I ||ki − k′i|| << min i,j∈I i 6=j ||ki − kj ||. We can calculate an under-approximation Approx∩(K,K ′) of the inter- section of two close polytopes K ∩K ′ using the following method: 1. For each corresponding pair of vertices (ki, k ′ i): (a) check whether ki ∈ K ′. If so, let k̃i = ki. (b) check whether k′i ∈ K. If so, let k̃i = k′i. (c) If neither 1a nor 1b apply, we set k̃i to be the point of K ∩ K ′ closest to the midpoint of ki and k ′ i. This point can be expressed as the solution of the following convex problem k̃i = arg min ||k − ki + k ′ i 2 || subject to k ∈ K ∩K ′. 2. Let Approx∩(K,K ′) = Conv({k̃i|i ∈ I}). Under-approximating the viability kernel Using the functions Approx⊕ and Approx∩, we can now give the following algorithm. This algorithm computes a polytopic under-approximation of the viability kernel which contains no more vertices than the viability constraint set. 34 Algorithm 3.2 Approximate polytopic method K0 ← K n← 0 while n ≤ N do if Kn = ∅ then KN ← ∅ break end if if Kn = Kn−1 then KN ← Kn break end if L← A−1 Approx⊕(Kn,V) . Approximation of Reach(Kn). Kn+1 ← Approx∩(Kn, L) . Approximation of Kn ∩ L. n← n+ 1 end while return (KN ) . KN ⊆ V iabT(K) 3.2.3 Ellipsoidal method Using ellipsoids provides an alternative approach to keeping the complex- ity of our set representation constant. This leads to a scalable method of computing an under-approximation of the viability kernel over large time horizons in high-dimensional spaces. The issue with this algorithm is that the class of ellipsoids is closed under neither Minkowski sum nor intersec- tion. Hence both must be under- approximated, leading to a reduction in accuracy. The Ellipsoidal Toolbox [23] provides a set of algorithms for under- approximating the set reachable from an ellipsoidal initial set under ellip- soidal constraints and computing an ellipsiodal under-approximation of the intersection of two ellipsoids. Details of these algorithms are provided in [19]. We use ApproxReach(K, `) to denote the under-approximation of the backward reachable set over a single discrete time step Reach1(K) tight in 35 the direction `, and use Approxmaxvol(K,L) to denote the maximum-volume ellipsoid contained in the intersection of K and L. Using this notation, we present Algorithm 3.3 which computes an under-approximation of the via- bility kernel using an ellipsoidal set representation. Algorithm 3.3 Ellipsoidal method K0 ← K n← 0 while n ≤ N do if Kn = ∅ then KN ← ∅ break end if if Kn = Kn−1 then KN ← Kn break end if L← ApproxReach(Kn, `) . Under-approximation of Reach(Kn). Kn+1 ← Approxmaxvol(Kn, L) . Under-approximation of Kn ∩ L. n← n+ 1 end while return (KN ) . KN ⊆ V iabT(K) 3.2.4 Support vector method We now present a method of under-approximating the viability kernel using support vectors. This method is based on ideas developed in [25] which uses support functions to compute an over-approximation of reachable sets. The computation of the viability kernel presents an additional challenge compared with the computation of reach sets in that intersections must also be performed. We first present a method of computing the support function of the viability kernel in a given direction by finding the solution to a convex 36 optimization problem. The solution to this optimization can then be used to compute a support vector in the same direction by means of a recursive formula. Approximating the viability kernel using support functions and support vectors provides an advantage over ellipsoids in terms of accuracy. Given an arbitrary direction ` ∈ Rd, the support function method allows us to find a hyperplane tangent to the viability kernel in the direction `. Performing this computation in multiple directions `k allows us to over-approximate the (discrete-time) viability kernel as the intersection of half-spaces bounded by the tangent hyperplanes. Similarly, we can compute a tight under- approximation of the viability kernel as the convex hull of a set of support vectors in the directions `k. This procedure can be made arbitrarily accurate simply by choosing a sufficient number of directions `k. The support vector method has an advantage over the polytope method in terms of scalability. After presenting Algorithm 3.4, we demonstrate its scalability experimentally using a chain of integrators of varying length. By Theorem 6 we can express the value of the support function ofKn+1 = A−1(Kn ⊕ V) ∩K0 in the direction ` as ρKn+1(`) = inf w∈Rd {ρK0(`− w) + ρV(A−Tw) + ρKn(A −Tw) } (3.7) where A−T = (AT )−1 is the inverse transpose of A. The function w 7→ ρK0(` − w) + ρV(A−Tw) + ρKn(A−Tw) is convex so if the functions ρK0 , ρV and ρKn could be evaluated in constant time, this problem could be solved efficiently. However, since ρKn in turn depends on ρKn−1 , a naive implementation of this formula would result in a number of calls to ρK0 that is exponential in n. This problem can be avoided by writing a closed-form expression for ρKn . Theorem 8. The value of the support function ρKn in the direction ` can be expressed as the solution to a convex optimization over an nd-dimensional 37 space. It is given by ρKn(`) = inf w∈Rnd ξ(`, w) where ξ(`, w) = ξ(`, w1, · · · , wn) = ρK0(`− wn) + n−1∑ k=1 ρK0(A −Twk+1 − wk) (3.8) +ρK0(A −Tw1) + n∑ k=1 ρV(A−Twk). Proof. Follows from (3.7) by induction on n. Using Theorem 7 we can express the set of support vectors of the set Kn+1 in the direction ` as vKn+1(`) = vA−1(Kn⊕V)∩K0(`) (3.9) = vK0(`− w̄) ∩A−1 ( vV(A−T w̄)⊕ vKn(A−T w̄) ) where w̄ ∈ arg inf w∈Rd { ρK0(`− w) + ρV(A−Tw) + ρKn(A−Tw) } . We have the following algorithm for computing under-approximations P↓ and over-approximations P↑ of the viability kernel that are tight in the set of directions L. The set P↓ is an under-approximation of the true via- bility kernel, hence all states in P↓ are guaranteed to be viable under the constraints K. The over-approximation (which may contain states that are not viable under K) is used to provide an upper bound on the error in the under-approximation. 38 Algorithm 3.4 Support vector method for ` ∈ L do minimize ξ(`, w) = ρK0(`− wn) + n−1∑ k=1 ρK0(A −Twk+1 − wk) +ρK0(A −Tw1) + n∑ k=1 ρV(A−Twk). subject to w ∈ Rnd ρ(`)← ξ(`, w̄) . Minimum value stored as ρ(`) V0 = vK0(A −T w̄n) for k = 1 . . . n do Vk = vK0(A −T w̄n−k − w̄n−k+1) ∩ A−1 ( vV(A−T w̄n−k+1)⊕ Vk−1 ) end for v(`)← Vn end for P↓ ← Conv (⋃ `∈L v(`) ) P↑ ← ⋂ `∈L {x | x · ` ≤ ρ(`)} return (P↓, P↑) . P↓ ⊆ V iabn(K0) ⊆ P↑ The following proposition simplifies the computation of a support vector vKn(`) in a number of practical cases. For example, it is used in the imple- mentation given in the supplementary materials for the case when the input and state constraints are ellipsoidal. Theorem 9. Suppose that Kn is nonempty. Let w̄ ∈ arg minw ξ(`, w) where 39 ξ is defined as in Theorem 8. If for some natural number k < n the set A−kvK0(A −T w̄k − w̄k+1)⊕ k⊕ i=1 A−ivV(A−T w̄i) is a singleton then vKn(`) = A −kvK0(A −T w̄k − w̄k+1)⊕ k⊕ i=1 A−ivV(A−T w̄i). Proof. Using (3.10, we can express vKn(`) as vKn(`) = vK0(`− w̄1) ∩A−1 ( V(A−T w̄1)⊕ ( vK0(A −T w̄1 − w̄2) ∩A−1 ( vV(A−T w̄2)⊕ · · · ⊕ vK0(A−T w̄n−1 − w̄n) ∩A−1 ( vV(A−T w̄n)⊕ vK0(A−T w̄n) ) · · · ))) ⊆ A−kvK0(A−T w̄k − w̄k+1)⊕ k⊕ i=1 A−ivV(A−T w̄i). However sinceKn is compact and nonempty we know that vKn(`) is nonempty. Thus the fact that A−kvK0(A −T w̄k − w̄k+1)⊕ k⊕ i=1 A−ivV(A−T w̄i) is a singleton implies that equality must hold. 3.3 Comparison of algorithms We now compare the four Lagrangian algorithms that we have presented. The exact polytope, approximate polytope, ellipsoid and support vector algorithms are implemented in MATLAB using the Multi-Parametric Tool- box v. 2.6.3 [24] for the polytope-based Algorithms 3.1 and 3.2, Ellipsoidal Toolbox v. 1.1.3 [23] for the ellipsoidal Algorithm 3.3 and CVX [15] for the support vector-based Algoritm 3.4. Level set approximations to the viability kernel are computed using the Level Set Toolbox v. 1.1 [31]. Computations 40 were performed using MATLAB release 2009b on a machine with an Intel Pentium 4 processor running at 3.00 GHz and 2GB RAM. The MATLAB code to generate the figures in this chapter and Chapter 4 can be downloaded from the web at http://www.ece.ubc.ca/∼jmaidens/thesis matlab.zip 3.3.1 Accuracy Algorithm 3.1 computes the viability kernel exactly but Algorithms 3.2–3.4 rely on under-approximations either to reduce their computational burden, or because of the limitations of the chosen set representation. The approxi- mation accuracy is affected by the number of time steps (hence the number of approximations that must be made) and can also be improved by adjust- ing parameters like the number of vertices of the polytope representation, the number of ellipsoids in the ellipsoidal representation and the number of directions |L| in the support vector representation. We compare the accuracy of the three approximate algorithms by com- paring their performance on a standard example: the double integrator. Consider the discrete-time double integrator [ x1(t+ 1) x2(t+ 1) ] = [ 1 ρ 0 1 ][ x1(t) x2(t) ] + [ 1 2ρ 2 ρ ] u(t) u(t) ∈ U = [−u0, u0]. (3.10) Figures 3.3 – 3.5 compare the accuracy of the three approximate algorithms as the time horizon is varied using the model (3.10) with a time step of ρ = 0.1 time units and input constraint |u(t)| ≤ 0.3. Figures 3.6 – 3.8 compare the accuracy of the three approximate algorithms as the number of vertices, ellipsoids and support vectors are varied using a time step of ρ = 0.1 time units, a horizon of 20 time steps and input constraint |u(t)| ≤ 0.3. 41 −0.4 −0.2 0 0.2 0.4 0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 x1 x 2 5 time steps −0.4 −0.2 0 0.2 0.4 0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 x1 x 2 10 time steps −0.4 −0.2 0 0.2 0.4 0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 x1 x 2 15 time steps −0.4 −0.2 0 0.2 0.4 0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 x1 x 2 30 time steps Figure 3.3: Viability constraints K = {x : ||x||∞ ≤ 0.5} (light blue) and the corresponding under-approximation of the viability ker- nel V iab(K) (grey) for the double integrator (3.10) computed using Algorithm 3.2. The approximation is performed using a fixed-complexity polytope with 8 vertices. 42 −0.4 −0.2 0 0.2 0.4 0.6 −0.4 −0.2 0 0.2 0.4 0.6 5 time steps −0.4 −0.2 0 0.2 0.4 0.6 −0.4 −0.2 0 0.2 0.4 0.6 10 time steps −0.4 −0.2 0 0.2 0.4 0.6 −0.4 −0.2 0 0.2 0.4 0.6 15 time steps −0.4 −0.2 0 0.2 0.4 0.6 −0.4 −0.2 0 0.2 0.4 0.6 30 time steps Figure 3.4: Viability constraints K = {x : ||x||2 ≤ 0.5} (light blue) and a corresponding single-ellipsoid under-approximation of the viability kernel V iab(K) (grey) for the double integrator (3.10) computed using Algorithm 3.3. The level set approximation (which is a highly accurate approximation of the continuous- time viability kernel, and is computationally feasible for this low-dimensional system) is shown in dark blue. 43 −0.4 −0.2 0 0.2 0.4 0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 x1 x 2 5 time steps −0.4 −0.2 0 0.2 0.4 0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 x1 x 2 10 time steps −0.4 −0.2 0 0.2 0.4 0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 x1 x 2 15 time steps −0.4 −0.2 0 0.2 0.4 0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 x1 x 2 30 time steps Figure 3.5: Viability constraints K = {x : ||x||2 ≤ 0.5} (light blue) and a corresponding support vector under-approximation of the viability kernel V iab(K) (grey) for the double integrator (3.10) computed using Algorithm 3.4. The level set approximation (which is a highly accurate approximation of the continuous- time viability kernel, and is computationally feasible for this low-dimensional system) is shown in dark blue. 44 −0.4 −0.2 0 0.2 0.4 0.6 −0.4 −0.2 0 0.2 0.4 0.6 x1 x 2 Polytope with 4 vertices −0.4 −0.2 0 0.2 0.4 0.6 −0.4 −0.2 0 0.2 0.4 0.6 x1 x 2 Polytope with 8 vertices −0.4 −0.2 0 0.2 0.4 0.6 −0.4 −0.2 0 0.2 0.4 0.6 x1 x 2 Polytope with 16 vertices −0.4 −0.2 0 0.2 0.4 0.6 −0.4 −0.2 0 0.2 0.4 0.6 x1 x 2 Polytope with 32 vertices Figure 3.6: Viability constraints K = {x : ||x||∞ ≤ 0.5} (light blue) and the corresponding under-approximation of the viability ker- nel V iab(K) (grey) for the double integrator (3.10) computed using Algorithm 3.2. The approximation is performed using a polytopes with varying numbers of vertices. 45 −0.4 −0.2 0 0.2 0.4 0.6 −0.4 −0.2 0 0.2 0.4 0.6 Ellipsoids computed in 4 directions −0.4 −0.2 0 0.2 0.4 0.6 −0.4 −0.2 0 0.2 0.4 0.6 Ellipsoids computed in 8 directions −0.4 −0.2 0 0.2 0.4 0.6 −0.4 −0.2 0 0.2 0.4 0.6 Ellipsoids computed in 16 directions −0.4 −0.2 0 0.2 0.4 0.6 −0.4 −0.2 0 0.2 0.4 0.6 Ellipsoids computed in 32 directions Figure 3.7: Viability constraints K = {x : ||x||2 ≤ 0.5} (light blue) and a corresponding piecewise-ellipsoidal under-approximation of the viability kernel V iab(K) (grey) for the double integra- tor (3.10) computed using Algorithm 3.3. The level set ap- proximation (which is a highly accurate approximation of the continuous-time viability kernel, and is computationally feasible for this low-dimensional system) is shown in dark blue. 46 −0.4 −0.2 0 0.2 0.4 0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 x1 x 2 Support vector computed in 4 directions −0.4 −0.2 0 0.2 0.4 0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 x1 x 2 Support vector computed in 8 directions −0.4 −0.2 0 0.2 0.4 0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 x1 x 2 Support vector computed in 16 directions −0.4 −0.2 0 0.2 0.4 0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 x1 x 2 Support vector computed in 32 directions Figure 3.8: Viability constraints K = {x : ||x||2 ≤ 0.5} (light blue) and a corresponding support vector under-approximation of the viability kernel V iab(K) (grey) for the double integrator (3.10) computed using Algorithm 3.4. The level set approximation (which is a highly accurate approximation of the continuous- time viability kernel, and is computationally feasible for this low-dimensional system) is shown in dark blue. 47 In Figures 3.3–3.5 we see that the approximate polytope algorithm (Al- gorithm 3.2) loses accuracy fastest as the number of time steps increases. The piecewise-ellipsoidal method (Algorithm 3.3) also loses accuracy as the number of time steps increases, but to a lesser extent. The accuracy of the support vector method (Algorithm 3.4) does not depend on the number of time steps. In Figures 3.6–3.8 we see that the quality of the approximation can be improved by increasing the size of the set representation. The quality of the polytope approximation is improved significantly by increasing the number of vertices from 4 to 8, but the improvement beyond 8 vertices is negligible. The accuracy of the piecewise-ellipsoidal algorithm continues to improve as the number of ellipsoids increases, but is still limited by the fact that all approximations must be ellipsoid-shaped. The support vector algorithm has the advantage that it can be made arbitrarily precise by choosing a sufficient number of directions in which to compute a support vector. Overall, the support vector algorithm performs best in terms of accu- racy, followed by the ellipsoidal algorithm, with the approximate polytopic method performing worst. In the next section, we compare the performances of the two most accurate approximation-based algorithms, along with the polytope method given in Algorithm 3.1 which computes the viability kernel exactly. 3.3.2 Scalability We compare how well the three algorithms scale as a function of the state dimension by comparing their performance on the chain of n integrators ẋ(t) = 0 1 0 · · · 0 0 0 1 0 ... . . . 0 0 0 1 0 0 0 · · · 0 x(t) + 0 0 0 1 u(t) u(t) ∈ U = [−u0, u0] (3.11) 48 discretized with a time step of ρ = 0.4. Its viability kernel is computed over a horizon of 10 steps with input constraint |u(t)| ≤ 0.3 and state constraint set K = {x : ||x||∞ ≤ 0.5} for the polytope algorithm and K = {x : ||x||2 ≤ 0.5} for the ellipsoid and support vector algorithms. In Figure 3.9 we plot the time it takes to compute (a) the viability kernel using the polytope method (b) an ellipsoidal under-approximation to the viability kernel in a single direction (c) 2n support vectors on the boundary of the viability kernel tight in a set of directions consisting of the standard basis vectors in Rn and their negatives. 5 10 15 20 25 30 35 400 20 40 60 80 100 120 140 160 180 200 state dimension run tim e ( s) polytope ellipsoid support vector Figure 3.9: Comparison of the run time for a chain of integrators of length n. 3.4 High-dimensional example The forced heat equation ∂ξ ∂t (x, t) = α ∂2ξ ∂x2 (x, t)− β(ξ − ξ0) + u(x, t) 49 describes how the temperature distribution over a finite, one-dimensional rod evolves over time when a heat input u is provided to the system. We assume that heat is lost from the rod at a rate (with proportionality constant β) dependent on the difference between the rod’s temperature and the ambient temperature ξ0 (assumed to be zero), that the system’s thermal diffusivity is α and that the rod is heated from one end. To study solutions to this equation numerically, we consider the temperature ξi at d discrete points xi on a uniform one-dimensional lattice of spacing ∆x. After finite difference approximation of the spatial derivative, we get the dynamics ξ̇(t) = − ( α 2(∆x)2 L+D ) ξ(t) +Bu(t) where ξ is now the vector of temperatures at the lattice points, L is the lattice’s Laplacian matrix, D is the diagonal matrix with entries β and B = [1 0 0 · · · 0]T . We compute the viability kernel for this system over a time horizon τ = 5 for the state constraint set consisting of a sphere of radius 10 centred at [10, . . . , 10]T using a lattice of 20 points, yielding a 20-dimensional system. We set α = 8, β = 0.015, ∆x = 1 and constrain the input u to the set [0, 5]. We compute an approximation of the viability kernel in a set of 120 directions consisting of 12 uniformly-spaced unit vectors in the ξ2k−1 × ξ2k plane for k = 1, . . . , 10. Figure 3.10 shows the resulting approximation of the viability kernel projected onto a selection of coordinate planes. Note that under- and over-approximations appear tight in the projec- tions onto the ξ1 × ξ2 and ξ19 × ξ20 planes since we sampled more support vectors in these subspaces. The projections onto the ξ1 × ξ10 and ξ1 × ξ20 subspaces could be improved simply by sampling in more directions. 3.5 Summary In this chapter, we presented four Lagrangian algorithms for approximat- ing viability kernels. These algorithms allow us to analyse systems with larger state dimension than was previously feasible using existing Eulerian 50 0 5 10 15 200 5 10 15 20 ξ1 ξ 2 0 5 10 15 200 5 10 15 20 ξ19 ξ 20 0 5 10 15 200 5 10 15 20 ξ1 ξ 10 0 5 10 15 200 5 10 15 20 ξ1 ξ 20 Figure 3.10: Viability kernel for the heat conduction problem pro- jected onto various coordinate planes. The projection of the state constraint set is shown in red and under- and over- approximations of the viability kernel are shown in blue and grey respectively. This approximation in 120 directions was performed in t = 425.16 seconds. As expected, we see that the set of viable states appears smaller when projected onto coordinates farther from the heated end at ξ1. techniques. The polytopic method given in Algorithm 3.1 computes the vi- ability kernel exactly, but quickly becomes infeasible as the state dimension increases. Conversely, the approximate polytopic method given in Algorithm 3.2 scales well with the state dimension, but exhibits limited accuracy. The ellipsoidal method of Algorithm 3.3 strikes a balance between scalability and accuracy, allowing an under-approximation of the viability kernel to be computed for high-dimensional systems with relatively good accuracy. The 51 support vector method given in Algorithm 3.4 performs the best in terms of both scalability and accuracy. It scales better than the previous algorithms (As demonstrated in Figure 3.9), and it can be made arbitrarily accurate as the number of support vector evaluations increases. 52 Chapter 4 Fallback mode initiation for physiologic closed-loop controllers Closed-loop controllers are becoming increasingly commonplace in modern healthcare. They have the potential to make medical treatment more pre- dictable, safer and less costly. However, they also pose potential risks to patient safety if the automation does not function as intended. Thus the regulatory approval process is expected to be the greatest challenge in get- ting these new devices to market [29]. The International Electrotechnical Commission (IEC) has published a collateral standard, IEC 60601-1-10:2007, to address issues related to phys- iologic closed-loop controllers [17]. Among other things, this standard spec- ifies recommendations for fallback modes and new verification/validation requirements for physiologic closed-loop control algorithms beyond the nor- mal verification standards for medical device software. In this chapter, we consider the problem of determining when the closed-loop device should initiate its fallback mode. Fallback mode is intended to provide a backup mode of operation should the closed-loop control system fail to work as intended, for example if it allows a physiologic variable to leave the range designated as safe. We 53 approach the fallback mode initiation problem from a viability theory per- spective. This method allows us to predict in advance when a physiologic variable will inevitably leave the safe range despite any admissible control input that could be applied, allowing the fallback mode to be initiated before any physiologic variable leaves the safe range. 4.1 Constraints for physiologic variables 4.1.1 Input and state constraints The behaviour of a physiologic closed-loop control system (PCLCS) is gov- erned by the dynamic relationship between a manipulated variable (or in- put) u that is applied to the system by a physiologic closed-loop controller (PCLC) and a set of measured or estimated physiologic variables (or states) x that are to be controlled in order to match a set of reference variables. To ensure safe operation of the PCLCS, the physiologic variables are de- sired to remain within a safe operating range (or viability constraint set) K. This range can be determined experimentally and could ideally be personal- ized to individual patients based on demographic and medical information. If the state of the PCLCS leaves K, we consider a fault to have occurred, causing an alarm to sound and fallback mode to be initiated. Additional constraints U are placed on the input variable, either because of safety-related concerns or due to limitations of the actuators that manip- ulate the system’s input. These constraints are not considered useful for fault detection and fallback mode initiation but rather describe limitations on the controllability of the PCLCS. 4.1.2 Refining constraints for early fallback mode initiation Impending violations of the viability constraints can be detected early by considering the dynamic relationships among the state and input variables. If it is determined that no admissible input will be able to keep the physio- logic variables within the constraints, a fallback can be initiated before the viability constraints are violated. This ability for early fault detection will 54 improve the safe operations of physiologic closed-loop control devices. The dynamic relationship between input and state variables under con- straints is captured by the viability kernel. The viability kernel is a refined subset of the set of viability constraints that takes the input constraints and system dynamics into account. Its significance to the fallback mode initi- ation problem is illustrated in the next section, using a problem from the closed-loop control of mechanical ventilation. 4.2 Illustrative example: closed-loop mechanical ventilation We consider the two-dimensional model of the end-tidal carbon dioxide par- tial pressure response to minute ventilation described in [16]: ẋ1 = k12(x2 − x1)− k10x1 − 1V1u ẋ2 = k21(x1 − x2) (4.1) where z1 = x1 + b represents the alveolar CO2 partial pressure (in the lungs) and z2 = x2 + b represents the venous CO2 partial pressure (in the body tissues). Carbon dioxide diffuses between the compartments with rate constants k12 and k21 and is eliminated from the body through the lungs at a rate dependent on k10 and the minute ventilation input u. The transfer diagram for this model is shown in Figure 4.1. The empirically-determined model parameters, determined based on the nominal values given in [16], are given in Table 4.1. Figure 4.1: Transfer diagram for the respiratory system model (4.1). 55 Table 4.1: Indentified parameters for the model (4.1) Parameter Identified value V1 45.4 L k10 1.08× 10−2 min−1 k21 2.61× 10−2 min−1 k12 4.93× 10−2 min−1 b 10.2 kPa The manipulated variable u (minute ventilation) is constrained by its nature to be positive and we specify for reasons of safety that the ventilation rate not exceed a maximum value umax. We wish that the venous partial pressure z2 remain within a narrow range c ≤ z2 ≤ C. Since there are no explicit safety concerns related to the alveolar CO2 partial pressure, we do not put an explicit constraint on the first state z1, other than that it be nonnegative. However, due to the diffusion of blood gasses between the two compart- ments, the constraints on z2 induce implicit constraints on z1. To make this argument concrete, we will study numerical solutions to (4.1) under the constraints c = 2 kPa and C = 10 kPa. In Figure 4.2 we see that even under the maximal input u(t) = umax = 10 L/min initializing the system with z1 = 20 kPa and z2 = 9 kPa leads to the violation of the constraint z2 ≤ C. Thus the point z = (20, 9) lies outside the viability kernel. That is, even though the constraints 2 ≤ z2 ≤ 10 are not violated by the initial point z = (20, 9), we should initiate the fallback mode in advance because the state will eventually leave the constraint set regardless of any admissible input that we apply. When the fallback mode is initiated, the ventilator automatically switches to an open-loop mode of operation and the attend- ing health professional is immediately notified that the device is not longer operating in closed-loop mode. The right-most boundary of the viability kernel is provided by the tra- jectory that attains a maximum z2 value at z2 = C under the maximal 56 Figure 4.2: Phase plane for system (4.1) with u(t) = umax. The tra- jectory forming the right-most boundary of the viability kernel is shown in bold. control input. This trajectory is plotted in bold in the phase plane diagram shows in Figure 4.2. The bottom-right-most corner of the constraint set is also not viable because the blood gasses in the two compartments equilibrate quicker than gas can build up due to natural processes. Even when no respi- ratory input is applied (u(t) = 0), a trajectory of (4.1) with initial condition z1 = 0 kPa and z2 = 2.1 kPa will violate the constraint z2 ≥ c because the venous partial pressure z2 will dip below c = 2 as it equilibrates with the arterial partial pressure. This left-most boundary trajectory is shown in bold in Figure 4.3. The viability kernel for this system is computed using the polytope method described in Algorithm 3.1 and shown in Figure 4.4. 57 Figure 4.3: Phase plane for system (4.1) with u(t) = 0. The trajectory forming the lower-left-most boundary of the viability kernel is shown in bold. 58 0 5 10 15 20 25 30 35 40 2 3 4 5 6 7 8 9 10 z1 z 2 Figure 4.4: Viability kernel for system (4.1) computed using Algo- rithm 3.1. 59 4.3 Closed-loop control of anaesthesia The Electrical & Computer Engineering in Medicine group at the Univer- sity of British Columbia is currently conducting preliminary clinical tests of a paediatric closed-loop anaesthesia system at the British Columbia Chil- dren’s Hospital. To ensure safe operation, this system places hard bounds on the Propofol administration rate (the input) and compartmental Propofol concentrations (the states). In this context, we compute the viability kernel for the purposes of “fall- back mode” initiation [17]. If at any point in time, the control system determines that the system’s state is outside the viability kernel, an alarm should be sounded and a fallback mode initiated. When the fallback mode is initiated, the Propofol infusion is immediately halted to mitigate any risk of overdose. The attending anaesthetist is provided information regarding the fault that has occurred (e.g. whether the Propofol concentration is too high or too low, and in which compartment the violation has/will occur) and he is instructed to begin a manual titration. We use the three-compartment model ċ1(t)ċ2(t) ċ3(t) = −(k10 + k12 + k13) k12 k13k21 −k21 0 k31 0 −k31 c1(t)c2(t) c3(t) + 1/V10 0 u(t− td) u(t) ∈ U = [0, u0] with input delay td = 0.5 minutes and the model parameters for an 11 year-old child of 35 kg taken from the Paedfusor data set [1]. The delay is approximated using a third order Padé approximation and then the system is discretized with a time step ρ = 0.25 minutes, yielding a 6-dimensional discrete-time model. We simulate a 90 minute surgery (time horizon τ = 90 minutes) with input constraint u(t) ∈ [0, 7 000] µg/min (200µg/kg/min) to compute the viability kernel of the state constraints given by an ellipsoid 60 (a) Viability constraints (red) and under-approximation of the viability kernel (blue) (b) Viability constraints (red) and over-approximation of the vi- ability kernel (grey) Figure 4.5: Viability kernel for the 6-dimensional pharmacokinetic model given in Subsection 4.3 computed using Algorithm 3.4. Constraints and viability kernels are shown projected onto the first 3 coordinates of the 6-dimensional state space. These ap- proximations were computed in t = 693s. inscribed in the box [1, 6] × [0, 10] × [0, 10] × [−100, 100] × [−100, 100] × [−100, 100]. The computed approximation is sampled in the set of directions L = Q −1 ` 0 0 0 : ` ∈ VI ∪ Q −1 0 0 0 ±ei : i = 1, . . . , 3 where VI is the set of the 12 vertices of an icosahedron in R3 and ei are the standard basis vectors in R3. The computed over- and under-approximations of the viability kernel are shown in Figure 4.5. We see that initial states with a low concentration in the slowly-equilibrating compartment c3 are not viable under the delayed dynamics. 61 4.4 Summary In this chapter we discussed the problem of determining when a closed-loop physiologic control system should initiate a fallback mode of operation, and we presented a viability-theoretic solution to this problem. We illustrated our approach with two applications: a closed-loop mechanical ventilation system and the closed-loop administration of the anaesthetic Propofol. 62 Chapter 5 Conclusions In this thesis, we presented a novel connection between viability and reach- ability that enables us to compute viability kernels in terms of backward reachable sets. This allowed us to develop and compare four algorithms that approximate the viability kernel using Lagrangian methods. The algo- rithm based on support vectors performed exceedingly well in the sense that it is scalable as the number of time steps or the state dimension increases and it can be made arbitrarily accurate. We also argued the applicability of viability theory for the advanced detection of safety violations in physiologic closed-loop control systems. 5.1 Future research 5.1.1 Differential games Throughout this thesis, we considered the problem of computing viability kernels, which correspond to the case where all inputs to the system are co- operative, or controllable. In the case of competing inputs and disturbances, the analogous construction is known as the discriminating kernel, defined as: Disc(K) = {x ∈ K : ∃u ∈ U ∀d ∈ D x(0) = x⇒ ∀t ≥ 0 x(t) ∈ K} The extension of some of the algorithms given in Section 3.2 to the com- 63 putation of discriminating kernels and reachable sets of differential games is straight-forward, while other algorithms (the support function method in particular) have no straight-forward extension. Extending the support vec- tor method would allow the computation of reach sets and discriminating kernels for high-dimensional linear differential games with convex constraint sets. 5.1.2 Improved computations with support vectors The support vector method described in Algorithm 3.4 relies on computing the intersection of numerous sets of support vectors. However, the issue of how to best represent, compute and intersect these sets numerically has not been discussed. Our current implementation relies on the shortcut given in Theorem 9, but a robust implementation of Algorithm 3.4 that applies to arbitrary linear systems with arbitrary convex constraint sets will require better numerical techniques for computing and representing sets of support vectors. 5.1.3 Time-scale separation methods Many high-dimensional system models, particularly those that arise in the study of biological systems, exhibit separations in time scale. The trajecto- ries of these systems converge quickly to an invariant submanifold (or linear subspace) of the state space. Singular perturbation methods can be used to decouple slow and fast time scales [18], allowing reachability or viability computations for a high- dimensional system to be performed in complimentary lower-dimensional spaces. If a high-dimensional system can be decomposed into systems of sufficiently small dimension, it would allow the use of a wider range of tech- niques (e.g. level set methods) for the computation of higher-dimension viability kernels. 64 5.1.4 Application to anaesthesia The application of our fallback mode initiation criterion discussed in Chap- ter 4, while mathematically solid, remains a theoretical approach. It still remains to implement a robust version of Algorithm 3.4 and integrate it with the closed-loop anaesthesia system developed by UBC’s Electrical & Computer Engineering in Medicine group. This would allow us to test the efficacy of our fallback mode initiation system in a clinical setting. 65 Bibliography [1] A. Absalom and G. Kenny. Paedfusor pharmacokinetic data set. British Journal of Anaesthesia, 95:110, 2005. → pages 60 [2] J.-P. Aubin. Viability Theory. Systems and Control: Foundations and Applications. Birkhäuser, Boston, MA, 1991. → pages 4, 10 [3] J.-P. Aubin and H. Frankowska. Set-Valued Analysis. Systems and Control: Foundations and Applications. Birkhäuser, Boston, MA, 1990. → pages 10 [4] J.-P. Aubin, A. M. Bayen, and P. Saint-Pierre. Viability Theory: New Directions. Springer, 2011. → pages 4, 9 [5] D. Avis and K. Fukuda. A pivoting algorithm for convex hulls and vertex enumeration of arrangements and polyhedra. Discrete & Computational Geometry, 8:295–313, 1992. → pages 27 [6] F. Blanchini and S. Miani. Set-Theoretic Methods in Control. Springer, 2008. → pages 1, 4, 32 [7] N. Bonneuil. Computing the viability kernel in large state dimension. Journal of Mathematical Analysis and Applications, 323:1444 – 1454, 2006. → pages 6 [8] L. Chapel, G. Deffuant, S. Martin, and C. Mullon. Defining yield policies in a viability approach. Ecological Modelling, 212:10–15, 2008. → pages 6 [9] A. Chutinan and B. Krogh. Computational techniques for hybrid system verification. Automatic Control, IEEE Transactions on, 48:64 – 75, 2003. → pages 1, 7, 25 [10] E. M. Clarke, O. Grumberg, and D. A. Peled. Model Checking. MIT Press, 2000. → pages 1 66 [11] P. Coquelin, S. Martin, and R. Munos. A dynamic programming approach to viability problems. In IEEE International Symposium on Approximate Dynamic Programming and Reinforcement Learning, pages 178 –184, 2007. → pages 6 [12] G. Deffuant, L. Chapel, and S. Martin. Approximating viability kernels with support vector machines. IEEE Transactions on Automatic Control, 52(5):933–937, 2007. → pages 5 [13] G. Frehse, C. Le Guernic, A. Donzé, S. Cotton, R. Ray, O. Lebeltel, R. Ripado, A. Girard, T. Dang, and O. Maler. SpaceEx: Scalable verification of hybrid systems. In Proceedings of the Conference on Computer Aided Verification, 2011. → pages 9 [14] A. Girard, C. Le Guernic, and O. Maler. Efficient computation of reachable sets of linear time-invariant systems with inputs. In J. Hespanha and A. Tiwari, editors, Hybrid Systems: Computation and Control, LNCS 3927, pages 257–271. Springer-Verlag, 2006. → pages 25 [15] M. Grant and S. Boyd. Graph implementations for nonsmooth convex programs. In Recent Advances in Learning and Control, pages 95–110. Springer-Verlag, 2008. → pages 40 [16] J.-O. Hahn, P. Shrawane, G. Dumont, and J. Ansermino. An empirical model of end-tidal CO2 response to minute ventilation. In Engineering in Medicine and Biology Conference, pages 524 –527, 2010. → pages 55 [17] IEC 60601-1-10:2007. Collateral standard: Requirements for the development of physiologic closed-loop controllers, 2007. → pages 53, 60 [18] C. K. R. T. Jones. A geometric approach to systems with multiple time scales. Bulletin of JSIAM (Ouyou-Suri), 7, 1997. → pages 64 [19] S. Kaynama, J. Maidens, I. M. Mitchell, M. Oishi, and G. A. Dumont. Computing the viability kernel using maximal reachable sets. In Hybrid Systems: Computation and Control, pages 55–64, 2012. → pages 10, 35 [20] A. B. Kurzhanski and P. Varaiya. Ellipsoidal techniques for reachability analysis. In N. Lynch and B. Krogh, editors, Hybrid 67 Systems: Computation and Control, LNCS 1790, pages 202–214, Berlin Heidelberg, 2000. Springer-Verlag. → pages 7, 9 [21] A. B. Kurzhanski and P. Varaiya. Ellipsoidal techniques for reachability analysis: internal approximation. Systems & Control Letters, 41:201–211, 2000. → pages 25 [22] A. B. Kurzhanski and P. Varaiya. Dynamic optimization for reachability problems. Journal of Optimization Theory and Applications, 108:227–251, 2001. → pages 9 [23] A. A. Kurzhanskiy and P. Varaiya. Ellipsoidal toolbox (ET). In Proceedings of IEEE Conference on Decision and Control, pages 1498–1503, San Diego, CA, Dec. 2006. → pages 28, 35, 40 [24] M. Kvasnica, P. Grieder, M. Baotić, and M. Morari. Multi-Parametric Toolbox (MPT). In R. Alur and G. J. Pappas, editors, Hybrid Systems: Computation and Control, LNCS 2993, pages 448–462, Berlin, Germany, 2004. Springer. → pages 28, 32, 40 [25] C. Le Guernic. Reachability analysis of hybrid systems with linear continuous dynamics. PhD thesis, Université Grenoble 1 - Joseph Fourier, 2009. → pages 4, 7, 29, 36 [26] C. Le Guernic and A. Girard. Reachability analysis of linear systems using support functions. Nonlinear Analysis: Hybrid Systems, 4(2): 250 – 262, 2010. ISSN 1751-570X. doi:10.1016/j.nahs.2009.03.002. IFAC World Congress 2008. → pages 25 [27] J. Lygeros. On reachability and minimum cost optimal control. Automatica, 40(6):917–927, June 2004. → pages 9 [28] J. Lygeros, C. J. Tomlin, and S. Sastry. Controllers for reachability specifications for hybrid systems. Automatica, 35:349–370, 1999. → pages 1 [29] P. J. Manberg, C. M. Vozella, and S. D. Kelley. Regulatory challenges facing closed-loop anesthetic drug infusion devices. Clin Pharmacol Ther, 84:166–169, 2008. → pages 53 [30] I. M. Mitchell. Comparing forward and backward reachability as tools for safety analysis. In A. Bemporad, A. Bicchi, and G. Buttazzo, editors, Hybrid Systems: Computation and Control, LNCS 4416, pages 428–443, Berlin Heidelberg, 2007. Springer-Verlag. → pages 13 68 [31] I. M. Mitchell and J. A. Templeton. A toolbox of Hamilton-Jacobi solvers for analysis of nondeterministic continuous and hybrid systems. In M. Morari and L. Thiele, editors, Hybrid Systems: Computation and Control, LNCS 3414, pages 480–494, Berlin, Germany, 2005. Springer-Verlag. → pages 9, 40 [32] I. M. Mitchell, A. M. Bayen, and C. J. Tomlin. A time-dependent Hamilton-Jacobi formulation of reachable sets for continuous dynamic games. IEEE Transactions on Automatic Control, 50(7):947–957, July 2005. → pages 5 [33] D. Panagou, K. Margellos, S. Summers, J. Lygeros, and K. J. Kyriakopoulos. A viability approach for the stabilization of an underactuated underwater vehicle in the presence of current disturbances. In Proceedings of IEEE Conference on Decision and Control, pages 8612–8617, Dec. 2009. → pages 4 [34] P. A. Parrilo. Structured Semidefinite Programs and Semialgebraic Geometry Methods in Robustness and Optimization. PhD thesis, California Institute of Technology, 2000. → pages 4 [35] S. Prajna and A. Jadbabaie. Safety verification of hybrid systems using barrier certificates. In Hybrid Systems: Computation and Control, pages 477–492. Springer, 2004. → pages 3 [36] R. M. Rifkin and R. A. Lippert. Value regularization and fenchel duality. Journal of Machine Learning Research, 8:441–479, 2007. → pages 72 [37] R. Rockafellar and R. Wets. Variational analysis. Grundlehren der mathematischen Wissenschaften. Springer, 1998. → pages 29, 71, 72 [38] P. Saint-Pierre. Approximation of the viability kernel. Applied Mathematics and Optimization, 29(2):187–209, Mar 1994. → pages 5, 14, 16, 17 [39] C. J. Tomlin, I. M. Mitchell, A. M. Bayen, and M. Oishi. Computational techniques for the verification and control of hybrid systems. Proceedings of the IEEE, 91(7):986–1001, 2003. → pages 1, 4, 9 69 Appendix A Proof of Theorem 7 As the proof of Theorem 7 requires the introduction of some elementary convex analysis, we have placed it here in the appendix. A.1 Convex conjugation and duality For a proper convex function f : Rd → R̄ = R∪{∞}, its conjugate function f∗ : Rd → R̄ is a proper convex function defined as f∗(`) = sup x∈Rd {〈x, `〉 − f(x)}. If f is further assumed to be lower semi continuous (lsc) then we have f∗∗ = f and hence there is a duality, or conjugacy correspondence, between lsc proper convex functions. Define the characteristic function δS : Rd → R̄ of a nonempty convex set S ⊂ Rd δS(x) = { 0 if x ∈ S ∞ if x 6∈ S The convex conjugate of the characteristic function δS is the support func- 70 tion σS σS(`) = sup x∈S 〈x, `〉 = sup x∈Rd {〈x, `〉 − δS(x)} = δ∗S(`). The conjugation operation also induces a duality between operations on functions. Define the infimal convolution (or epi-sum) f1#f2 of two convex functions f1 and f2 as f1#f2 (`) = inf w∈Rd f1(`− w) + f2(w). The following proposition establishes that infimal convolution is the dual of the summation operation. Lemma 1 (Theorem 11.23a, Rockafellar and Wets [37]) If f1 and f2 are two lsc proper convex functions whose domains have nonempty intersection then • (f1 + f2)∗ = f∗1 # f∗2 • (f1 # f2)∗ = f∗1 + f∗2 . A.2 Subgradients and support vectors Let f : Rd → R̄ = (−∞,∞] be convex. The set ∂f(`) of subgradients of f at ` is the set of vectors x ∈ Rd such that for all ˜̀∈ Rd f(`)− f(˜̀) ≥ x · (`− ˜̀) The following result is standard: Lemma 2 (Fenchel’s Inequality, Proposition 11.3, Rockafellar and Wets [37]) Let f : Rd → R̄ be a proper convex function. Then for all `, x ∈ Rd we have f(`)− ` · x+ f∗(x) ≥ 0. Further, equality holds if and only if x ∈ ∂f(`). From Fenchel’s Inequality, we get the following characterization of the set of support vectors of a compact, convex set U in the direction `: 71 Theorem 10. Let σU (`) and vU (`) be defined as in (3.3) and (3.5). Then vU (`) = ∂σU (`) Proof. x ∈ vU (`) ⇐⇒ x ∈ U and x · ` = σU (`) ⇐⇒ σU (`)− x · `+ δU (x) = 0 ⇐⇒ x ∈ ∂σU (`). A.3 Theorem 7 Before proving Theorem 7 we need one final result: Lemma 3 (Proposition 2.22a and Theorem 10.13, Rockafellar and Wets [37]) Let g : Rd × Rn → R̄ be a proper convex function. Then g̃(`) = infw∈Rn g(`, w) is a proper convex function and (x, 0) ∈ ∂g(`, w̄)⇐⇒ x ∈ ∂g̃(`) and g̃(`) = g(`, w̄). The first two bullet points in Theorem 7 follow directly from the defini- tion. We prove only the third. The proof is adapted from [36]. Define g(`, w) = σU (` − w) + σV(w) so that g̃(`) = infw∈Rn g(`, w) = (σU#σV)(`). The assumption w̄ = arg infw∈Rd{σU (` − w) + σV(w)} means 72 that g̃(`) = g(`, w̄). Hence it follows from Lemmas 2 and 3 that x ∈ vU∩V(`) = ∂σU∩V(`) = ∂(σU#σV)(`) = ∂g̃(`) ⇐⇒ (x, 0) ∈ ∂g(`, w̄) ⇐⇒ 0 = g(`, w̄)− x · `− 0 · w̄ + g∗(x, 0) = σU (`− w̄) + σV(w̄)− x · `+ σ∗U (x) + σ∗V(x+ 0) = [σU (`− w̄)− x · (`− w̄) + σ∗U (x)] +[σV(w̄)− x · w̄ + σ∗V(x)] ⇐⇒ x ∈ ∂σU (`− w̄) and x ∈ ∂σV(w̄) ⇐⇒ x ∈ vU (`− w̄) ∩ vV(w̄). 73
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Scalable computation of viability kernels and a viability-theoretic...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Scalable computation of viability kernels and a viability-theoretic approach to guaranteeing safety for… Maidens, John Norman 2012
pdf
Page Metadata
Item Metadata
Title | Scalable computation of viability kernels and a viability-theoretic approach to guaranteeing safety for closed-loop medical devices |
Creator |
Maidens, John Norman |
Publisher | University of British Columbia |
Date Issued | 2012 |
Description | As closed-loop controllers become increasingly prevalent in medical technology, increasing emphasis is being placed on ensuring that such systems operate in a safe manner. In our approach to guaranteeing the safe operation of a physiologic closed-loop control system, we wish to provide a mathematical guarantee that, despite limited control authority, the system’s state can be confined to a region designated as safe. The largest subset of the safe region for which there exists an admissible control input that keeps the state within the safe region is known as the viability kernel, or maximal controlled invariant set. Many methods are known for computing viability kernels in low-dimensional systems, but these existing methods rely on gridding the state space and hence their time complexity increases exponentially with the state dimension. In this thesis we describe a new connection between reachability and viability theory that enables us to approximate the viability kernel using Lagrangian methods which scale well with the state dimension. We present four new viability kernel approximation algorithms using polytope-, ellipsoid- and support vector-based set representations and we compare their performances in terms of accuracy and scalability with the state dimension. Using the support vector and ellipsoidal techniques, we are able to accurately approximate the viability kernel for systems of much larger state dimension than was previously feasible using existing Eulerian methods. We also present a viability theoretic solution to the problem of determining when a physiologic closed-loop control system should initiate a fallback mode of operation. The viability-based method allows impending safety violations to be detected in advance, allowing the fallback mode to be initiated earlier than using a naive approach. Our new approach to fallback mode initiation is examined in two sample contexts: the closed-loop control of carbon dioxide partial pressure under mechanical ventilation, and the control of the concentration of the anaesthetic drug Propofol using a paediatric model of Propofol pharmacokinetics. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2012-07-16 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-ShareAlike 3.0 Unported |
DOI | 10.14288/1.0072881 |
URI | http://hdl.handle.net/2429/42719 |
Degree |
Master of Applied Science - MASc |
Program |
Biomedical Engineering |
Affiliation |
Applied Science, Faculty of |
Degree Grantor | University of British Columbia |
GraduationDate | 2012-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-sa/3.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2012_fall_maidens_john.pdf [ 1.42MB ]
- Metadata
- JSON: 24-1.0072881.json
- JSON-LD: 24-1.0072881-ld.json
- RDF/XML (Pretty): 24-1.0072881-rdf.xml
- RDF/JSON: 24-1.0072881-rdf.json
- Turtle: 24-1.0072881-turtle.txt
- N-Triples: 24-1.0072881-rdf-ntriples.txt
- Original Record: 24-1.0072881-source.json
- Full Text
- 24-1.0072881-fulltext.txt
- Citation
- 24-1.0072881.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0072881/manifest