Scalable Techniques for the Computation of Viable and Reachable Sets Safety Guarantees for High-Dimensional Linear Time-Invariant Systems by Shahab Kaynama A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in The Faculty of Graduate Studies (Electrical and Computer Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) July 2012 c© Shahab Kaynama 2012 Abstract Reachability analysis and viability theory are key in providing guarantees of safety and proving the existence of safety-preserving controllers for con- strained dynamical systems. The minimal reachable tube and (by duality) the viability kernel are the only constructs that can be used for this purpose. Unfortunately, current numerical schemes that compute these constructs suffer from a complexity that is exponential in the dimension of the state, rendering them impractical for systems of dimension greater than three or four. In this thesis we propose two separate approaches that improve the scal- ability of the computation of the minimal reachable tube and the viability kernel for high-dimensional systems. The first approach is based on structure decomposition and aims to facilitate the use of computationally intensive yet versatile and powerful tools for higher-dimensional linear time-invariant (LTI) systems. Within the structure decomposition framework we present two techniques—Schur-based and Riccati-based decompositions—that im- pose an appropriate structure on the system which is then exploited for the computation of our desired constructs in lower-dimensional subspaces. The second approach is based on set-theoretic methods and draws a new connection between the viability kernel and maximal reachable sets. Exist- ing tools that compute the maximal reachable sets are efficient and scalable with polynomial complexity in time and space. As such, these scalable tech- niques can now be used to compute our desired constructs and therefore provide guarantees of safety for high-dimensional systems. Based on this new connection between the viability kernel and maximal reachable sets we propose a scalable algorithm using ellipsoidal techniques for reachability. We show that this algorithm can efficiently compute a conservative under- ii Abstract approximation of the viability kernel (or the discriminating kernel when uncertainties are present) for LTI systems. We then propose a permissive state-feedback control strategy that is capable of preserving safety despite bounded input authority and possibly unknown disturbances or model un- certainties for high-dimensional systems. We demonstrate the results of both of our approaches on a number of practical examples including a problem of safety in control of anesthesia and a problem of aerodynamic flight envelope protection. iii Preface Part of the work presented in this thesis has been published in or submitted to several international conferences and scientific journals. Results from Chapter 3 were published in: [i] S. Kaynama and M. Oishi, “Complexity reduction through a Schur- based decomposition for reachability analysis of linear time-invariant systems,” International Journal of Control, vol. 84, no. 1, pp. 165–179, January 2011 [ii] S. Kaynama and M. Oishi, “Schur-based decomposition for reacha- bility analysis of linear time-invariant systems,” in Proc. IEEE Con- ference on Decision and Control, Shanghai, China, December 16–18, 2009, pp. 69–74 Results from Chapter 4 will appear in: [iii] S. Kaynama and M. Oishi, “A modified Riccati transformation for complexity reduction in reachability analysis of linear time-invariant systems,” IEEE Transactions on Automatic Control (accepted as Full Paper) Results from Chapter 5 were published in: [iv] S. Kaynama, J. Maidens, M. Oishi, I. M. Mitchell, and G. A. Du- mont, “Computing the viability kernel using maximal reachable sets,” in Proc. Hybrid Systems: Computation and Control, Beijing, China, April 17–19, 2012, pp. 55–63 Appendix C and Section 5.3.2 are partially based on the results that ap- peared in: iv Preface [v] S. Kaynama, M. Oishi, I. M. Mitchell, and G. A. Dumont, “The con- tinual reachability set and its computation using maximal reachability techniques,” in Proc. Joint IEEE Conference on Decision and Con- trol, and European Control Conference, Orlando, FL, December 12– 15, 2011, pp. 6110–6115 [vi] S. Kaynama, M. Oishi, I. M. Mitchell, and G. A. Dumont, “Fixed- complexity piecewise ellipsoidal representation of the continual reacha- bility set based on ellipsoidal techniques,” to appear in Proc. American Control Conference, Montreal, QC, June 27–29, 2012 In all papers I was the lead investigator. I determined the problem to solve, conducted the theoretical research and verified the results via simula- tions. I was responsible for the bulk of the writing of the papers. My super- visors and co-authors Professors M. Oishi, I. M. Mitchell, and G. A. Dumont provided technical and editing feedback. The paper [iv] is a collaborative work with J. Maidens. I formulated the problem, formed the hypothesis of expressing the viability kernel in terms of maximal reachable sets, and presented the solution for the discrete-time case with help and advice from Professor Mitchell. Mr. Maidens offered the solution for the continuous- time case (Propositions 5.1 and 5.2), while I helped with the proofs and the presentation of the results. I developed the computational algorithms and verified the methods via simulations. Professors Oishi, Mitchell, and Dumont provided technical and editing feedback throughout the work. v Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Preface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Reachability Analysis and Viability Theory . . . . . . . . . . 3 1.3 Computational Techniques & the “Curse of Dimensionality” 4 1.4 The Goal and Contributions of the Thesis . . . . . . . . . . . 5 1.4.1 Complexity Reduction via Structure Decomposition . 6 1.4.2 Complexity Reduction via Set-Theoretic Methods . . 7 1.5 Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.5.1 Complexity Reduction for Safety Analysis . . . . . . 8 1.5.2 Safety-Preserving Control Synthesis . . . . . . . . . . 10 1.6 Organization of the Thesis . . . . . . . . . . . . . . . . . . . 12 1.7 Basic Notation . . . . . . . . . . . . . . . . . . . . . . . . . . 13 vi Table of Contents 2 Preliminaries and Problem Formulation . . . . . . . . . . . 15 2.1 Backward Constructs for Constrained Dynamical Systems . . 15 2.2 Significance of Minimal Reach Tube and Viability Kernel . . 19 2.3 Main Goal and Problem Formulation . . . . . . . . . . . . . 21 2.4 Approach I: Structure Decomposition . . . . . . . . . . . . . 22 2.4.1 Objective and Problem Statement . . . . . . . . . . . 22 2.4.2 Definitions and Preliminaries . . . . . . . . . . . . . . 22 2.5 Approach II: Set-Theoretic Methods . . . . . . . . . . . . . . 25 2.5.1 Objective and Problem Statement . . . . . . . . . . . 25 2.5.2 Definitions and Preliminaries . . . . . . . . . . . . . . 25 2.6 Robust Reachability Analysis: Competing Inputs . . . . . . 26 2.7 To Over- or Under-Approximate? . . . . . . . . . . . . . . . 29 3 Schur-Based Structure Decomposition . . . . . . . . . . . . 30 3.1 Disjoint Control Input . . . . . . . . . . . . . . . . . . . . . . 31 3.2 Non-Disjoint Control Input . . . . . . . . . . . . . . . . . . . 33 3.3 Reachability in Lower Dimensions . . . . . . . . . . . . . . . 34 3.3.1 Formulating an Upper-Bound on Growth of Z1[0,τ ] . . 36 3.4 Further Reduction of Complexity . . . . . . . . . . . . . . . 38 3.5 Extension to Switched Linear Systems . . . . . . . . . . . . . 40 3.6 Numerical Examples . . . . . . . . . . . . . . . . . . . . . . . 41 3.6.1 Arbitrary 3D System . . . . . . . . . . . . . . . . . . 42 3.6.2 4D Aircraft Dynamics . . . . . . . . . . . . . . . . . . 44 3.6.3 8D Distillation Column . . . . . . . . . . . . . . . . . 44 3.6.4 4D Unstable System (An Example for Section 3.4) . . 47 3.7 Summary and Further Discussions . . . . . . . . . . . . . . . 48 4 Riccati-Based Structure Decomposition . . . . . . . . . . . . 51 4.1 Disjoint Control Input . . . . . . . . . . . . . . . . . . . . . . 53 4.2 Non-Disjoint Control Input . . . . . . . . . . . . . . . . . . . 53 4.2.1 Transformation i . . . . . . . . . . . . . . . . . . . . . 54 4.2.2 Transformation ii . . . . . . . . . . . . . . . . . . . . 58 4.2.3 The Unidirectional Coupling Term (Choosing δ) . . . 60 vii Table of Contents 4.3 Recursive Decomposition . . . . . . . . . . . . . . . . . . . . 62 4.4 Reachability in Lower Dimensions . . . . . . . . . . . . . . . 63 4.4.1 Formulating an Upper-Bound on Growth of Z2[0,τ ] . . 65 4.4.2 Dimension vs. Magnitude of Uncertainty . . . . . . . 67 4.4.3 Conservatism Due to Projection . . . . . . . . . . . . 67 4.4.4 Implications of Riccati-Based Reachable Tube . . . . 68 4.5 Numerical Examples . . . . . . . . . . . . . . . . . . . . . . . 70 4.5.1 Arbitrary 4D System . . . . . . . . . . . . . . . . . . 70 4.5.2 4D Cart with Two Inverted Pendulums . . . . . . . . 72 4.5.3 Arbitrary 6D System . . . . . . . . . . . . . . . . . . 74 4.5.4 Comparison With Schur-Based Decomposition . . . . 75 4.5.5 The Decomposition and Maximal Reachability . . . . 77 4.6 Summary and Further Discussions . . . . . . . . . . . . . . . 79 5 Set-Theoretic Methods: Lagrangian Tools for Viability . . 81 5.1 The Viability Kernel in Terms of Maximal Reach Sets . . . . 82 5.1.1 Continuous-Time Systems . . . . . . . . . . . . . . . 82 5.1.2 Discrete-Time Systems . . . . . . . . . . . . . . . . . 87 5.2 Computational Algorithms . . . . . . . . . . . . . . . . . . . 88 5.2.1 A Piecewise Ellipsoidal Approach . . . . . . . . . . . 89 5.3 Practical Examples . . . . . . . . . . . . . . . . . . . . . . . 97 5.3.1 Flight Envelope Protection (Continuous-Time) . . . . 97 5.3.2 Safety in Anesthesia Automation (Discrete-Time) . . 100 5.4 Summary and Further Discussions . . . . . . . . . . . . . . . 103 6 Robustness and Safety-Preserving Control Synthesis . . . 106 6.1 Discriminating Kernel vs. Robust Maximal Reach Sets . . . 106 6.1.1 Under-Approximating the Discriminating Kernel . . . 107 6.2 Computational Algorithm & Control Synthesis . . . . . . . . 109 6.2.1 Piecewise Ellipsoidal Approximation . . . . . . . . . . 110 6.2.2 Safety-Preserving Feedback Policy . . . . . . . . . . . 112 6.2.3 Remarks and Practical Considerations . . . . . . . . 118 6.3 Numerical Examples . . . . . . . . . . . . . . . . . . . . . . . 122 viii Table of Contents 6.3.1 2D System Without Uncertainty . . . . . . . . . . . . 122 6.3.2 2D System With Uncertainty . . . . . . . . . . . . . . 123 6.3.3 3D Control of Anesthesia . . . . . . . . . . . . . . . . 130 6.3.4 6D Flight Envelope Protection . . . . . . . . . . . . . 131 6.4 Summary and Further Discussions . . . . . . . . . . . . . . . 139 7 Conclusions and Future Work . . . . . . . . . . . . . . . . . . 143 7.1 Summary of Contributions . . . . . . . . . . . . . . . . . . . 144 7.2 Future Research Directions . . . . . . . . . . . . . . . . . . . 145 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 Appendices A Supplementary Materials for Chapter 3 . . . . . . . . . . . 160 A.1 On Assumption 3.1 (and 4.2) . . . . . . . . . . . . . . . . . . 160 A.2 Proof of Proposition 3.4 . . . . . . . . . . . . . . . . . . . . . 160 A.3 Decomposed System Matrices for Examples 3.6.2 and 3.6.3 . 164 A.3.1 4D Aircraft Dynamics (Example 3.6.2) . . . . . . . . 164 A.3.2 8D Distillation Column (Example 3.6.3) . . . . . . . 165 B Supplementary Materials for Chapter 4 . . . . . . . . . . . 166 B.1 Proofs of Propositions 4.2 and 4.3 . . . . . . . . . . . . . . . 166 B.2 Condition Number of the Modified Riccati Transformation . 167 B.3 Decomposed System Matrices for Example 4.5.3 and Section 4.5.5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168 B.3.1 Arbitrary 6D System (Example 4.5.3) . . . . . . . . . 168 B.3.2 Maximal Reachability Example (Section 4.5.5) . . . . 169 C Other Backward Reachability Constructs . . . . . . . . . . 170 C.1 Definitions and Connections . . . . . . . . . . . . . . . . . . 170 ix List of Tables 6.1 Model parameters for the patient (11 years old, 35 kg) . . . . 130 x List of Figures 1.1 Inadequate design may lead to constraint violations . . . . . . 2 2.1 The maximal reachable set Reach]t(K,U) . . . . . . . . . . . . 16 2.2 The maximal reachable tube Reach][0,τ ](K,U) . . . . . . . . . 17 2.3 The minimal reachable tube Reach[[0,τ ](K,U) . . . . . . . . . 18 2.4 The (finite-horizon) viability kernel V iab[0,τ ](K,U) . . . . . . 19 3.1 Phase-plane of a planar system with a non-convex target . . . 39 3.2 Eigenvalue scenarios violating Proposition 3.4 . . . . . . . . . 40 3.3 The non-convex target set in the transformed coordinates . . 42 3.4 Schur-based vs. actual reach tube for Example 3.6.1. . . . . . 43 3.5 Schur-based vs. actual viability kernel for Example 3.6.2 . . . 45 3.6 The Schur-based viability kernel of Example 3.6.3 . . . . . . . 47 3.7 Over-approximation of the reach tube for Example 3.6.4 (via Proposition 3.4) . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.1 Unidirectional coupling ‖δF (Z(δ))‖ for Example 4.5.1 . . . . 62 4.2 Unidirectional coupling vs. time for randomly generated sys- tems of dimension n . . . . . . . . . . . . . . . . . . . . . . . 68 4.3 Riccati-based vs. actual reach tube for Example 4.5.1 . . . . . 71 4.4 Riccati-based vs. actual viability kernel for Example 4.5.2 . . 73 4.5 Riccati-based reach tube for Example 4.5.3 . . . . . . . . . . 75 4.6 Fraction of randomized tests for which a solution existed . . . 76 4.7 The unit disc under Riccati- and Schur-based transformations 77 4.8 The Riccati decomposition for maximal reachability . . . . . 79 5.1 Iterative under-approximation of V iab[0,τ ](K,U) . . . . . . . . 84 xi List of Figures 5.2 Distance d > 0 from the boundary set ∂K. . . . . . . . . . . . 86 5.3 A finer time interval partition results in better approximation 95 5.4 Convergence plot of the error as a function of |P | . . . . . . . 95 5.5 Under-approximation of V iab[0,2](K,U) for Example 5.3.1 . . 98 5.6 Under-approximation of V iab[0,2](K,U) for Example 5.3.1 . . 99 5.7 Under-approximation of V iab[0,90](K,U) for Example 5.3.2 . . 104 6.1 The graph of the hybrid safety-preserving controller . . . . . 113 6.2 βα as a function of φ for a given design parameter α ∈ [0, 1). 121 6.3 Closed-loop trajectories for Example 6.3.1 . . . . . . . . . . . 124 6.4 Safety-preserving feedback policy for Example 6.3.1 . . . . . . 125 6.5 Closed-loop trajectories for Example 6.3.1 . . . . . . . . . . . 126 6.6 Safety-preserving control with less switching frequency . . . . 127 6.7 Closed-loop trajectories for Example 6.3.2 . . . . . . . . . . . 128 6.8 Safety-preserving feedback policy for Example 6.3.2 . . . . . . 129 6.9 Disturbance signal for Example 6.3.2 . . . . . . . . . . . . . . 129 6.10 Trajectories without safety control (Example 6.3.3) . . . . . . 131 6.11 Trajectories with smooth safety control (Example 6.3.3) . . . 132 6.12 Smooth safety-preserving feedback policy (Example 6.3.3) . . 132 6.13 uperf(t) as a percentage of u(t) in qperf (Example 6.3.3) . . . . 133 6.14 Active mode of H using the modified policy (Example 6.3.3) 133 6.15 Internal input γ using the modified policy (Example 6.3.3) . . 134 6.16 Location of x(t) within the domains of H (Example 6.3.3) . . 134 6.17 Trajectories with aggressive safety control (Example 6.3.3) . . 135 6.18 Aggressive safety-preserving feedback policy (Example 6.3.3) 135 6.19 Active mode of the hybrid automaton H (Example 6.3.3) . . 136 6.20 Internal input γ (Example 6.3.3) . . . . . . . . . . . . . . . . 136 6.21 Location of x(t) within the domains of H (Example 6.3.3) . . 137 6.22 Closed-loop trajectories w/o safety control (Example 6.3.4) . 140 6.23 Closed-loop trajectories with safety control (Example 6.3.4) . 141 A.1 Illustrating the non-negativity of the Hamiltonian . . . . . . . 164 xii Acknowledgments If I have seen further, it is by standing on the shoulders of giants. —Sir Isaac Newton This thesis would not have been possible without the support of a num- ber of people. I consider myself extremely lucky to have worked under the incredible supervision and mentorship of Professors Meeko Oishi, Ian Mitchell, and Guy Dumont. Their expertise in the field, their readiness to provide consistent feedback, and their willingness to allow me the freedom to work however, wherever, and on whatever topic I find fascinating have been instrumental in my development as a researcher. I would not be where I am today had it not been for their continuous support in every which way I can think of. And for all that I am forever grateful. I would like to express my gratitude to the rest of my examination committee Professors Ryozo Nagamune, Philip Loewen, Tim Salcudean, Shahriar Mirabbasi, Zinovy Reichstein, and Dušan Stipanović (U. of Illi- nois at Urbana-Champaign) for the time and effort spent in considering my candidacy. I am especially indebted to Professor Nagamune for providing great feedback on my research throughout these past few years. I have had wonderful collaborations with John Maidens—collaborations that I hope continue beyond our time at UBC. I have no doubt that a bright future awaits him. I would also like to take this opportunity to thank Dr. Alex A. Kurzhanskiy (UC Berkeley) for always generously and diligently answering my many questions regarding Ellipsoidal Toolbox. I have had the pleasure of being part of a great research group. Nikolai, Neda, Ahmad, Pouyan, and Pouria have all made my stay at UBC an inter- xiii Acknowledgments estingly nonlinear experience. I thank Nikolai Matni (Caltech) in particular for all our intellectually stimulating conversations (both before and after his departure from UBC) and our mutual fascination by the beauty of the field of controls. He has been both a fantastic friend and a great resource to bounce ideas off of. The Control Systems Reading Group meetings were inspiring and informative: Special thanks to Professors Bhushan Gopaluni, Farrokh Sassani, Ryozo Nagamune, and Jin-Oh Hahn (U. of Alberta) and the rest of the team, particularly Ehsan, Daniel, Devin, Marius, and Klaske. Many scholars at UBC, UMIST, and EMU have had lasting influence on my understanding of controls and mathematics. I would like to particularly thank Professor Runyi Yu (EMU) whose in-depth knowledge of the field, genuine care for his students’ success, and consistent support of my academic endeavors throughout the years have helped shape what I am today. Last but not least, I thank my family—Mom, Dad, Kaveh (and Kimia), and Sina—for their never-ending and unconditional love, encouragement, and moral support. Thank you for teaching me the importance of educa- tion and the gratification of success. I am who I am because of you. My wonderful wife, Blerina, to whom I wholeheartedly dedicate this thesis, has been incredibly supportive of me and my work throughout the past almost decade—since our very first date at Café Blue Note! Had it not been for her selfless sacrifices and her faith in me, this work would simply not have come to existence. Thank you sweetheart for helping me up when days were dark and cheering me on when they were not. I also thank Deni for putting up with me over the past couple of years, and for being the most well-behaved and disciplined teenager I have known. This work has been financially supported by NSERC Discovery Grants #327387 and #298211, NSERC Collaborative Health Research Project #CH- RPJ-350866-08, and the Institute for Computing, Information and Cognitive Systems (ICICS). xiv To Blerina. xv Chapter 1 Introduction 1.1 Motivation Constrained dynamical systems have deservingly received a tremendous amount of attention among researchers in both academia and industry due to the presence of safety constraints and hard bounds that appear in many practical scenarios. Providing guarantees of constraint satisfaction and fa- cilitating appropriate synthesis of constraint-satisfying controllers therefore is highly desirable, particularly in safety-critical applications. One example of such an application with safety constraints and bounded input authority is the closed-loop control of anesthesia [26]; safety (state or output constraints) may be defined in terms of prespecified therapeutic bounds on plasma concentration of the anesthetics, and the input (drug infu- sion rate) is physically bounded by the actuator limitations. To ensure safety of the patient and obtain regulatory certificates, guarantees of safety of the system can play an important role. While we postpone further discussions on the safety aspects of the control of depth of anesthesia to Sections 5.3.2 and 6.3.3, we shall point out that naive design of a controller—even when all constraints have been accounted for—may result in violation of the safety constraints. For example, Figure 1.1 shows the response of a patient using a model predictive controller (MPC) that has been designed with both input and safety (output in this case) constraints in mind. Additional pre-design analysis such as formulation of a terminal constraint set for the receding horizon optimization problem and/or other tricks and techniques that en- force invariance, (strong) feasibility, and other properties of the closed-loop system are warranted. The reader is referred to e.g. [11, 61, 84] for a thor- ough treatment of such techniques as well as other application areas where 1 1.1. Motivation 1000 2000 3000 4000 50000 0.2 0.6 1 1.4 O u tp u t (P se u d o- o cc u p an cy ) Time (s) Safety bounds Setpoint Response (a) Response of the closed-loop system. 0 1000 2000 3000 4000 5000 0.04 0.08 0.12 In p u t (I n fu si on R at e) Time (s) (b) Synthesized control law. Figure 1.1: Control of drug plasma concentration in a high-responder pa- tient (#80) undergoing 1.5 hr surgery using an inadequately- designed MPC. The input (drug infusion rate) is bounded above and below (u ∈ [0, 0.8]). Clearly, the hard constraint on the out- put (y ∈ [0.1, 1]) indicating the safety requirement is violated for nearly 20 mins of the surgery; The receding horizon optimization becomes infeasible with respect to its initial condition, resulting in “softening” of the safety constraint. constraint satisfaction in MPC is desirable. Another example is the so-called flight or aircraft envelope protection problem [81, 83, 116, 117], where the safety constraints are defined as the aircraft’s aerodynamic envelope and consequently the flight management system must ensure that certain combinations of states are avoided at all times to prevent stalling or other undesirable behaviors. Other application domains in which safety must be maintained despite bounded inputs include aircraft autolanders [8], collision avoidance [25, 28, 46, 89, 96], automated highway systems [80], control of under-actuated un- derwater vehicles [98], stockout prevention of constrained storage systems in manufacturing processes [13], management of a marine renewable resource [9], and safety in semi-autonomous vehicles [119], to name a few. 2 1.2. Reachability Analysis and Viability Theory 1.2 Reachability Analysis and Viability Theory Reachability analysis and viability theory provide solid frameworks for con- trol synthesis and trajectory analysis of constrained dynamical systems in a set-valued fashion (cf. [5, 12, 22, 67, 69, 118]) and have been utilized in diverse applications ranging from those listed above to the control of un- certain oscillatory systems [23] and verification of temporal properties of a toggle circuit [44], etc. Reachability analysis identifies the set of states forward (backward) reach- able by a constrained dynamical system from a given initial (target) set of states. Reachability analysis distinguishes itself from what simulations can achieve in the sense that with simulations a single trajectory or execution of the system corresponding to a single initial state is computed at a time, whereas with reachability analysis all points belonging to all possible tra- jectories or executions are computed at once from all possible initial states. Properties of the system inferred from such set-valued computations are therefore universal and hold true for the entire set of initial states. Simi- larly, viability theory provides a set-valued insight into the behavior of the trajectories inside a given constraint set. For example, the viability kernel is the set of initial states for which there exists at least one trajectory or execution of the input-constrained system that respects the state constraint for all time. The notions of maximal and minimal reachability analysis were intro- duced in [86]. Their corresponding constructs differ in how the time vari- able and the bounded input are quantified. In formation of the maximal reachability construct, the input tries to steer as many states as possible to the target set. In formation of the minimal reachability construct, the trajectories reach the target set regardless of the input applied. Based on these differences, the maximal and minimal reachable sets and tubes (the set of states traversed by the trajectories over the time horizon [70, 86]) are formed. The objects generated by each of these constructs have unique attributes: The maximal reachability constructs can be used to synthesize input policies that steer the trajectories of the system to the target (or, 3 1.3. Computational Techniques & the “Curse of Dimensionality” to analyze how the trajectories behave under uncertainty and/or external disturbance). The minimal reachability constructs, on the other hand, can be used to synthesize inputs that keep the trajectories of the system away from the target set. The viability kernel and the minimal reachable tube are duals of one another [18, 79], while the invariance kernel (the set of states that remain in the constraint set for all possible inputs for all time) is the dual of the maximal reachable tube [5, 69, 79]. It is shown in [86] and [5] that the minimal reachable tube and the viabil- ity kernel are the only constructs that can be used to prove safety/viability of the system and to synthesize inputs (controllers) that preserve this safety. As such it is highly desirable to be able to compute these constructs for possibly high-dimensional safety-critical constrained systems for which guarantees of safety and prevention of constraint violation are crucially important. 1.3 Computational Techniques and the “Curse of Dimensionality” In this thesis we are concerned with backward constructs generated by reach- ability analysis and viability theory. That is, we seek to compute the set of initial states that is formed under the input-constrained system for a given target/terminal or constraint set of states. In general an exact computation of the reachability or viability constructs is extremely difficult if not impossi- ble. Instead, approximations of these sets are computed. Such computations have historically been subject to Bellman’s “curse of dimensionality”; their complexity increases rapidly with the dimension of the continuous states [4]. Algorithms that approximate the backward constructs can be divided into two main categories [86]: Eulerian methods (e.g., [18, 33, 89, 106]) are capable of computing the viability kernel and by duality, the minimal reachable tube. Although versatile in terms of ability to handle complex dynamics and constraints, these algorithms rely on gridding the state space and therefore their computational complexity increases exponentially with the dimension of the state, rendering them impractical for systems of di- 4 1.4. The Goal and Contributions of the Thesis mension higher than three or four. The second category of algorithms are Lagrangian methods (e.g., [30, 39, 40, 50, 70, 73, 78]) that follow trajectories and compute the maximal reachable sets and tube in a scalable and efficient manner. These algorithms take advantage of compact set representations (e.g., ellipsoids and zonotopes) and/or in general the convexity of all con- straints.1 Their computational complexity is therefore usually polynomial in time and space, making them suitable for application to high-dimensional systems. 1.4 The Goal and Contributions of the Thesis In many applications such as those involving safety-critical systems the abil- ity to compute the minimal reachable tube and the viability kernel can be paramount: Guarantees of safety and synthesis of safety-preserving con- trollers can only be obtained using these constructs. However, if the system is of even moderate dimension, the Eulerian methods that to this date have exclusively2 been used to compute these constructs cannot be employed. This argument is the cornerstone of our research. This thesis presents our efforts to address the scalability aspect of the curse of dimensionality to enable the computation of the minimal reachable tube and the viability kernel for higher-dimensional systems. We do so using two separate approaches: • Complexity reduction via structure decomposition, our first approach, 1The convexity requirement may be circumvented if a non-convex constraint can be represented by the union of appropriate closed and bounded convex sets. In such a case the computations are performed using a few convex initial sets, which in turn affects the computational time only linearly in the number of initial sets. In the general case regardless of the convexity issues, the resulting set computed by Lagrangian methods over the entire time horizon is a close approximation of the (possibly non-convex) maximal reachable tube. 2One exception is the method of polytopes [76], a Lagrangian method developed from within the MPC community, for discrete-time LTI systems that allows for the computa- tion of the controlled-invariant subset (viability kernel) when all constraints are polytopes. However, similarly to Eulerian methods, this technique suffers from an exponential com- plexity for most systems: The number of vertices of the polytope increases rapidly with each successive Minkowski sum. The convex-hull operation performed by the technique is also known to be computationally demanding [93]. 5 1.4. The Goal and Contributions of the Thesis aims to facilitate the use of powerful Eulerian methods on higher- dimensional continuous-time continuously-valued linear time-invariant (LTI) systems (and by extension, hybrid systems with LTI continuous dynamics). As such, many of the benefits that Eulerian methods offer (e.g., ability to handle arbitrarily-shaped possibly unbounded non- convex constraints, or the synthesis of safety-preserving control laws) can now be taken advantage of for relatively higher-dimensional LTI systems. • Complexity reduction via set-theoretic methods, our second approach, bridges the gap between maximal and minimal constructs and aims to facilitate the use of the scalable and efficient Lagrangian methods for the computation of the viability kernel for continuous- and discrete- time (possibly nonlinear) systems. Thanks to these results we propose an efficient and scalable piecewise ellipsoidal algorithm that not only enables the approximation of the viability kernel for high-dimensional LTI systems, but also facilities a scalable synthesis of safety-preserving controllers. 1.4.1 Complexity Reduction via Structure Decomposition Decomposing the system into lower dimensional subsystems is an approach (referred to as structure decomposition) that has been utilized in the past in e.g. [88, 91, 113] as a means to reduce the computational complexity in reachability analysis. While applicable to nonlinear systems, these tech- niques assume that the system itself presents a certain structure that can be exploited. In this thesis we propose a number of techniques, based on the structure decomposition framework, applicable to LTI systems of generic form. We assume no initial structures. Our techniques (under certain condi- tions) will impose the desired structure on the system which is then exploited for decomposition purposes and ultimately for reduction of complexity in reachability analysis. The decomposition allows for the computation of the reachable tube in lower dimensions, thus reducing the computational burden significantly. 6 1.4. The Goal and Contributions of the Thesis We accomplish this through transformation of the system into appropriate coordinate spaces in which reachability analysis could be performed in lower- dimensional subspaces and is guaranteed to yield an over-approximation of the actual minimal reachable tube in that space. By performing the analysis in lower dimensions we obtain significant reduction in the computational costs, albeit at the expense of over-approximation. While also applicable to maximal reachable tube computation, our pro- posed techniques are of primary benefit to the computation of the minimal reachable tube (or by duality, the viability kernel) for which Eulerian meth- ods have traditionally been used exclusively. 1.4.2 Complexity Reduction via Set-Theoretic Methods Through our second approach we will show that the viability kernel can be expressed as a nested sequence of maximal reachable sets. By bridging the gap between the viability kernel and the maximal reachable sets we pave the way for more efficient computation of the viability kernel through the use of Lagrangian algorithms. Significant reduction in the computational costs can be achieved since instead of a single calculation with exponential complexity one can perform a series of calculations with polynomial complexity. Based on the well-studied ellipsoidal techniques for maximal reachability analysis [70] we propose an algorithm that computes a guaranteed piecewise ellipsoidal under-approximation of the viability kernel for LTI systems. We show that the proposed algorithm is efficient and scalable and can be used to address the safety needs of high-dimensional safety-critical systems (such as the anesthesia automation or the flight envelope protection problems mentioned at the beginning of this chapter). The presented connection between the maximal reachable sets and the viability kernel facilitates a more scalable computation of the viability kernel (and by duality, the minimal reachable tube) as well as the respective safety- preserving control laws. 7 1.5. Related Work 1.5 Related Work Complexity reduction for reachability analysis has been addressed by a num- ber of researchers. Methods to compute the reachability constructs (max- imal or minimal) for higher-dimensional systems can be divided into three main categories. First are techniques that take advantage of certain repre- sentations of sets [30, 40, 64, 65, 74, 76, 78, 108]. Second are techniques that make use of model approximation [42, 47], hybridization [3], projec- tion [44, 91] and structure decomposition [48, 113, 121]. Finally, third are methods that combine the approaches from the first two categories; For ex- ample, the proposed technique in [49] employs Krylov subspace projection along with low-dimensional polytopes to compute the maximal reachable sets/tube for large-scale affine systems. Related works on complexity reduc- tion for the computation of the minimal reachable tube and the viability kernel are surveyed next. Efficient techniques for synthesis of maximal reachability controllers that steer the system to a given target set have been extensively studied in the past. As the main objective of this thesis is safety and safety-preserving control, we refrain from surveying these techniques and only provide a few references: [39, 40, 64, 68, 69, 94]. Instead, here we will discuss relevant techniques that are capable of directly synthesizing safety-preserving control laws. 1.5.1 Complexity Reduction for Computation of Minimal Reachable Tube and Viability Kernel A projection scheme in [91] based on Hamilton-Jacobi (HJ) partial differen- tial equations (PDEs) over-approximates the projection of the actual mini- mal reachable tube in lower dimensional subspaces, with the unmodeled di- mensions treated as a disturbance. Similarly, [113] decomposes a full-order nonlinear system into either disjoint or overlapping subsystems and solves multiple HJ PDEs in lower dimensions. The computed minimal reachable tube for each subsystem is an over-approximation of the projection of the full-order minimal reachable tube onto the subsystem’s subspace. 8 1.5. Related Work More recently, in [88] a mixed implicit-explicit HJ formulation of the minimal reachable tube is presented for nonlinear systems whose state vector contains states that are integrators of other states. It is shown that the computational complexity of this new formulation is linear in the number of integrator states, while still exponential in the dimension of the rest of the state space. In [21] an approximate dynamic programming technique is proposed that, although still grid-based, enables a more efficient computation of the viability kernel. The viability kernel (similarly to [79]) is expressed as the zero sub-level set of the value function of the corresponding optimal control problem. The authors assume that the value function, which is a viscos- ity solution of a HJB PDE, is differentiable everywhere on the constraint set. The HJB PDE is then discretized and the resulting value function is numerically computed on a grid using a function approximator such as the k-nearest neighbor algorithm. An error bound on the (over-)approximation of the viability kernel is provided. The approximation converges to the true kernel in the limit, as the number of grid points goes to infinity. In [25] it is shown that for systems with order-preserving dynamics when the control set is a direct product of two compact intervals in R (i.e. the con- trol set is a two-dimensional, axis-aligned rectangle), the minimal reachable tube equals the intersection of two reachable tubes each with a constant in- put. The constant inputs take value on the opposite vertices of the input set where the maximum of one interval meets the minimum of the other. The reachable tubes with constant inputs can be computed efficiently resulting in an efficient computation of the original minimal reachable tube. This scheme also yields an efficient synthesis of safety-preserving control laws. Another related approach is the search for a barrier certificate [99] (a Lyapunov-like function) for the system ẋ = f(x, u), x ∈ X , u ∈ U , that forms a separating surface between any two given sets A and B. If there exists a function B : X → R such that B(x) ≤ 0 ∀x ∈ A, B(x) > 0 ∀x ∈ B, and LfB(x) ≤ 0 ∀(x, u) ∈ X × U along the zero level set of B (here Lf denotes the Lie derivative along f), then no trajectories will ever go from A to B. This technique can be adapted to analytically formulate the boundary 9 1.5. Related Work of the infinite-horizon viability kernel as well—if it is non-empty. The Lie derivative on the surface of the candidate certificate must now satisfy ∃u ∈ U ∀x ∈ X LfB(x) ≤ 0. It is shown in [99] that for systems with polynomial vector fields and semi-algebraic constraints (e.g. polynomial inequalities), efficient techniques based on Sum of Squares and semi-definite programming can be used to find the barrier certificate.3 1.5.2 Safety-Preserving Control Synthesis Safety-preserving controllers, the control policies associated with viability kernels and minimal reachable tubes, are capable of keeping the trajectories of the system within the safe region of the state space. Their synthesis has therefore received significant attention among researchers as a means to guarantee the safety/viability of constrained dynamical systems [2, 5, 11, 81, 116] (cf. [7, 12, 89] for more recent expositions). A classical approach is to use the information of the shape of the com- puted set (viability kernel or minimal reachable tube) by choosing the safety- preserving optimal control laws based on the contingent cone [14] or the proximal normal [5, 18] at every point on the boundary of the set as per Nagumo’s theorem and its generalizations (cf. [5]). Similarly, a given viabil- ity kernel can be used as a feasible or terminal constraint set to guarantee (strong) feasibility of the receding horizon optimization problem in an MPC framework [11, 61, 84, 100, 101].4 This results in (recursive) constraint sat- isfaction of the closed-loop system and therefore the generated control laws can be regarded as safety-preserving. In both cases discussed above the com- putational complexity of synthesizing such control laws is at best equal to that of computing the corresponding viability kernel or minimal reachable 3This method cannot be used to formulate the finite-horizon viability kernel which may be useful when e.g. the infinite-horizon kernel is empty. Moreover, there are no guarantees that a barrier certificate can be found (even with simple, stable linear dynamics for which a Lyapunov function can be systematically constructed). 4While alternative techniques exist that eliminate the need for terminal constraint sets for feasibility (such as ensuring that the horizon of the receding horizon optimization problem is “sufficiently” large), in this thesis we are only concerned with computing the viability kernel and the minimal reachable tube as well as their corresponding safety controllers and will not cover such techniques. 10 1.5. Related Work tube. Another approach is related to the previously described barrier certifi- cate technique. The technique has been extended in [102] to synthesize safety-preserving controllers using density functions and convex optimiza- tion methods for polynomial nonlinear systems with semi-algebraic con- straints and simple magnitude bounded inputs. The pros and cons of this approach are similar to those of computing a barrier certificate with the exception that for the synthesis problem a gridding of the state space is also required. Unfortunately, this gridding renders the technique intractable in high dimensions. A classification technique based on support vector machines (SVMs) is presented in [24] that approximates the viability kernel and yields an analytical expression of its boundary. A sequential minimal optimization algorithm is solved to compute the SVM that in turn forms a barrier function [95] (not to be confused with a barrier certificate mentioned above) on or close to the boundary of the viability kernel in the discretized space. While the method successfully reduces the computational time for the synthesis of control laws when the dimension of the input space is high, its applicability to systems with high-dimensional state space is limited as it relies on the same state space gridding approach used in the classical techniques e.g. [106]. Furthermore, the method does not provide any guarantees that the synthesized control laws are safety-preserving. The notion of approximate bisimulation [41] can be used to construct a discrete abstraction of the continuous state space such that the observed be- havior of the corresponding abstract system is “close” to that of the original continuous system. Girard et al. in a series of papers [17, 37, 38] use this notion to construct safety-preserving controllers for approximately bisimi- lar discrete abstractions and prove that the controller is correct-by-design for the original systems. The technique is applied to incrementally stable switched systems (for which approximately bisimilar discrete abstractions of arbitrary precision can be constructed) with autonomous or affine dynamics, and safety-preserving switched controllers are synthesized. The abstraction, however, relies on sampling of time and space (i.e. gridding) and therefore 11 1.6. Organization of the Thesis its applicability appears to be limited to low-dimensional systems. 1.6 Organization of the Thesis The remainder of the thesis is organized as follows. Chapter 2 provides necessary preliminaries and formally formulates the objectives of the the- sis. Chapter 3 describes our first technique, Schur-based decomposition, within the structure decomposition framework. Our second decomposition technique, Riccati-based decomposition, is presented in Chapter 4 building upon the results of Chapter 3 and providing a more in-depth treatment of the structure decomposition framework for complexity reduction in reach- ability analysis. These two chapters are based on the assumption that the system dynamics are LTI and aim to facilitate the applicability of compu- tationally intensive Eulerian methods for higher-dimensional LTI systems. Chapter 5 describes a set-theoretic approach that enables the use of effi- cient Lagrangian methods for the computation of the viability kernel, cir- cumventing entirely the need to employ computationally intensive Eulerian methods. An efficient algorithm is presented that approximates the viabil- ity kernel of LTI systems in a scalable fashion. Chapter 6 extends these results to systems with disturbances and then formulates a scalable and robust safety-preserving feedback control strategy for LTI systems. Chap- ter 7 summarizes the thesis reiterating its major contributions, and provides directions for future research. Finally, supplementary materials are provided in the Appendices: Ap- pendix A contains the proof of Proposition 3.4 (page 38) and elaborates on Assumptions 3.1 (page 34) and 4.2 (page 54) used in Chapters 3 and 4. Appendix B contains the proofs of Propositions 4.2 (page 60) and 4.3 (page 61) and also provides an upper-bound on the condition number of the proposed modified Riccati transformation matrix in Chapter 4. Appendix C is complementary to the reachability constructs used within the body of the thesis and presents additional backward constructs that are formed under constrained dynamical systems. Their connections to one another and to the constructs used within the thesis are formalized, aiming to help the reader 12 1.7. Basic Notation attain a more complete picture of the contributions of this thesis. 1.7 Basic Notation The quantifiers ∃ and ∀ are existential and universal, respectively. We denote the N×N identity matrix by IN and the inner product by 〈·, ·〉. For brevity, ‖·‖ denotes the infinity norm. For a constant matrix A = [aij ] ∈ Rm×n the induced norm is ‖A‖ := sup v∈Rn, v 6=0 ‖Av‖ ‖v‖ = max1≤j≤n m∑ i=1 |aij |. (1.1) For a Lebesgue measurable function f : R → Rn defined over an interval [ta, tb] we denote ‖f‖ := ‖f(·)‖L∞[ta,tb] = sup t∈[ta,tb] ‖f(t)‖ <∞. (1.2) The Hausdorff distance between any two nonempty subsets X and Y of a metric space (Rn, d) is defined as distH(X ,Y) := max { sup x∈X inf y∈Y d(x, y), sup y∈Y inf x∈X d(x, y) } . (1.3) The erosion of X by Y (also known as the Pontryagin difference between X and Y) is defined as X Y := {x | x+ y ∈ X ∀y ∈ Y}. (1.4) The Minkowski sum of X and Y is defined as X ⊕ Y := {x+ y | x ∈ X , y ∈ Y}. (1.5) The projection of a set A ⊆ X × Y onto X is defined as ProjX (A) := {x ∈ X | ∃y ∈ Y s.t. (x, y) ∈ A}. (1.6) 13 1.7. Basic Notation We denote by ◦X the interior of X , by ∂X its boundary, and by clX := ◦X∪∂X its closure. The set B(δ) denotes a norm-ball of radius δ ∈ R+ about the origin in Rn. The set X c denotes the complement of X in Rn. 14 Chapter 2 Preliminaries and Problem Formulation 2.1 Backward Constructs for Constrained Dynamical Systems Consider a continuously-valued dynamical system L(x(t)) = f(x(t), u(t)), x(0) = x0 (2.1) with state space X := Rn (a finite-dimensional vector space), state vector x(t) ∈ X , and input u(t) ∈ U where U is a compact and convex subset of Rm. Depending on whether the system evolves in continuous time (t ∈ R+) or discrete time (t ∈ Z+), L(·) denotes the derivative operator or the unit forward shift operator, respectively. The vector field f : X × U → X is assumed to be continuous in its arguments in the discrete-time case, and Lipschitz in x and continuous in u in the continuous-time case. Let U[0,t] := {u : [0, t]→ Rm measurable, u(t) ∈ U a.e.}. (2.2) With an arbitrary, finite time horizon τ > 0, for every t ∈ [0, τ ], x0 ∈ X , and u(·) ∈ U[0,t], there exists a unique trajectory xux0 : [0, t]→ X that satisfies the initial condition xux0(0) = x0 and the differential/difference equation (2.1). When clear from the context, we shall drop the subscript and superscript from the trajectory notation. For a nonempty state constraint/target set K ⊂ X this thesis is primarily concerned with the following backward constructs, their implications, and 15 2.1. Backward Constructs for Constrained Dynamical Systems K Figure 2.1: The maximal reachable set Reach]t(K,U) (blue/plain) for the target set K (green/patterned). A few sample initial conditions and trajectories are also shown. their connections to one another. (For other backward constructs please see Appendix C.) Definition 2.1 (Maximal Reachable Set). The maximal reachable set at time t is the set of initial states for which there exists an input such that the trajectories emanating from those states reach K exactly at time t (Fig- ure 2.1): Reach]t(K,U) := {x0 ∈ X | ∃u(·) ∈ U[0,t], xux0(t) ∈ K}. (2.3) Definition 2.2 (Maximal Reachable Tube). The maximal reachable tube1 over the horizon [0, τ ] is the set of initial states for which there exists an input such that the trajectories emanating from those states reach K at some time t ∈ [0, τ ] (Figure 2.2): Reach][0,τ ](K,U) := {x0 ∈ X | ∃u(·) ∈ U[0,τ ], ∃t ∈ [0, τ ], xux0(t) ∈ K}. (2.4) 1Also known as the possible victory domain [18], the attainability tube [70], or the capture basin [6]. 16 2.1. Backward Constructs for Constrained Dynamical Systems K Figure 2.2: The maximal reachable tube Reach][0,τ ](K,U) (also includes K itself). A few sample initial conditions and trajectories are also shown. Definition 2.3 (Minimal Reachable Tube). The minimal reachable tube2 over the horizon [0, τ ] is the set of initial states for which for every input there exits a time t ∈ [0, τ ] such that the trajectories emanating from those states reach K at t (Figure 2.3): Reach[[0,τ ](K,U) := {x0 ∈ X | ∀u(·) ∈ U[0,τ ], ∃t ∈ [0, τ ], xux0(t) ∈ K}. (2.5) Definition 2.4 (Viability Kernel). The (finite-horizon) viability kernel3 of K is the set of all initial states in K for which there exists an input such that the trajectories emanating from those states remain in K for all time t ∈ [0, τ ] (Figure 2.4): V iab[0,τ ](K,U) := {x0 ∈ X | ∃u(·) ∈ U[0,τ ], ∀t ∈ [0, τ ], xux0(t) ∈ K}. (2.6) What differentiates the above constructs from one another is the type 2Also known as the certain victory domain [18], or the capture set in the differential games literature [81] (not to be confused with the capture basin). 3The infinite-horizon viability kernel V iabR+(K,U) is also known as the maximal controlled-invariant subset [11]. 17 2.1. Backward Constructs for Constrained Dynamical Systems K Figure 2.3: The minimal reachable tube Reach[[0,τ ](K,U) (also includes K itself). A few sample initial conditions and trajectories are also shown. and order of quantifiers operating on the time and input variables. These seemingly subtle differences generate fundamentally distinct sets (with prop- erties that are unique to each of them). In particular, the following inclusions hold. Proposition 2.1. V iab[0,τ ](K,U) ⊆ K ⊆ Reach[[0,τ ](K,U) ⊆ Reach][0,τ ](K,U). (2.7) Proof. That V iab[0,τ ](K) ⊆ K is well-known [5]. To showK ⊆ Reach[[0,τ ](K,U) take x0 ∈ K and let τ = 0. Thus xux0(0) = x0 for any u(·) ∈ U[0,τ ] which implies x0 ∈ Reach[[0,τ ](K,U). To prove Reach[[0,τ ](K,U) ⊆ Reach][0,τ ](K,U), take x0 ∈ Reach[[0,τ ](K,U). We have that ∀u(·) ∈ U[0,τ ] ∃t ∈ [0, τ ] xux0(t) ∈ K =⇒ ∃u(·) ∈ U[0,τ ] ∃t ∈ [0, τ ] xux0(t) ∈ K ⇐⇒ x0 ∈ Reach][0,τ ](K,U). 18 2.2. Significance of Minimal Reach Tube and Viability Kernel K Figure 2.4: The (finite-horizon) viability kernel V iab[0,τ ](K,U) (blue/plain). A few sample initial conditions and trajectories are also shown. 2.2 Significance of Minimal Reachable Tube and Viability Kernel What is so significant about computing the minimal reachable tube and the viability kernel? The short answer is “guarantees of safety” and “safety- preserving control synthesis”. To see this let us begin with maximal reach- ability constructs. The maximal reachable tube and sets can be used to synthesize input policies that steer the trajectories of the system to the target, or if the target is deemed “unsafe”, they can be used to analyze the set of states backward reachable in the worst case by the system under a bounded disturbance or uncertainty. The recently developed Lagrangian methods (e.g. [30, 39, 40, 50, 70, 73]) approximate the maximal reachable tube by first fixing the time variable and computing the corresponding maximal reachable set, and then taking the 19 2.2. Significance of Minimal Reach Tube and Viability Kernel union of these sets over the time horizon. It is easy to verify that Reach][0,τ ](K,U) ≡ ⋃ t∈[0,τ ] { x0 ∈ X | ∃u(·) ∈ U[0,t], xux0(t) ∈ K } (2.8) = ⋃ t∈[0,τ ] Reach]t(K,U). (2.9) The equality holds since both quantifiers operating on the input and time variables in (2.4) are existential and therefore their order can be inter- changed. This connection, which has been recognized and extensively used since the earlier works in the literature (e.g. [20, 22, 69]), was also formal- ized in [86] and proven using the Hamilton-Jacobi-Bellman framework in [79]. The Lagrangian methods are scalable and computationally efficient (with polynomial complexity in both time and space). On the other hand, the minimal reachable tube and the viability kernel can be used to synthesize safety-preserving inputs that keep the trajectories of the system away from the unsafe target set, or contained within a safe constraint set. In fact, it is shown in [86] and [5] that the minimal reachable tube, and by duality the viability kernel, are the only constructs that can be employed to prove the existence of an input which guarantees safety of the system. Notice that, as shown in [18] and as a dual of the results in [79], we have V iab[0,τ ](Kc,U) = (Reach[[0,τ ](K,U))c. (2.10) Since the viability kernel and the minimal reachable tube are duals of one another, they need not be treated separately. In this thesis we will initially focus on computing the minimal reachable tube for an unsafe target K in Chapters 3 and 4, and then shift focus to the case in which constraint K is deemed safe and therefore compute the viability kernel in the remaining chapters. Note however that approximating the minimal reachable tube and the viability kernel are not mutually exclusive—a method that facilitates an under-approximation of the viability kernel of K automatically provides a 20 2.3. Main Goal and Problem Formulation means for the over-approximation of the minimal reachable tube for Kc.4 In computing the minimal reachable tube the input is universally quan- tified. Furthermore, the time variable must be quantified only after the input is quantified. Interchanging the order of quantifiers here will change the meaning of the set and is therefore not possible. In fact, as also shown in [86], Reach[[0,τ ](K,U) ⊇ { x0 ∈ X | ∃t ∈ [0, τ ], ∀u(·) ∈ U[0,t], xux0(t) ∈ K } = ⋃ t∈[0,τ ] { x0 ∈ X | ∀u(·) ∈ U[0,t], xux0(t) ∈ K } . (2.11) Among Lagrangian methods, the technique in [72] has been extended to han- dle universally quantified inputs. However, since the time variable is quan- tified first, according to (2.11) the resulting set is an under-approximation of the unsafe minimal reachable tube. The powerful Eulerian methods (e.g. [18, 33, 89, 106]), in addition to computing the maximal reachable set, are capable of directly computing the minimal reachable tube and the viability kernel. However, they rely on gridding the state space and are therefore computationally intensive. Although versatile in terms of ability to handle various types of dynamics and constraints, the applicability of these techniques has been historically limited to systems of low dimensionality (up to 3D or 4D in practice) due to their exponential complexity. 2.3 Main Goal and Problem Formulation Given the fact that for guarantees of safety and the synthesis of safety- preserving controllers the minimal reachable tube and the viability kernel can play an extremely important role, we seek to devise efficient techniques that enable a more scalable computation of these constructs. To this end, 4Under-approximation of the viability kernel and over-approximation of the minimal reachable tube are the correct forms of approximation; See Section 2.7. 21 2.4. Approach I: Structure Decomposition we present two separate approaches. The first approach, the structure de- composition techniques that will be presented in Chapters 3 and 4, aims to enable the use of powerful Eulerian methods on higher-dimensional LTI sys- tems. The second approach, which is based on set-theoretic techniques, aims to facilitate the application of efficient and scalable Lagrangian methods for the computation of the viability kernel and safety-preserving controllers by drawing a connection between minimal and maximal reachability constructs. These objectives are defined next. 2.4 Approach I: Structure Decomposition 2.4.1 Objective and Problem Statement Consider the case in which (2.1) is a continuous-time LTI system ẋ = Ax+Bu (2.12) described by matrix notation S := [ A B ] (2.13) with A ∈ Rn×n and B ∈ Rn×p, and the constrained input u(t) ∈ U is the control input. Problem 2.1 (Objective I). Find an appropriate basis transformation for (2.12) such that in the new coordinates the system can be decomposed into lower-dimensional, decoupled or weakly-coupled subsystems for which reach- ability analysis can be performed independently and thus more efficiently. 2.4.2 Definitions and Preliminaries A linear transformation of the system S in (2.13) using a nonsingular matrix T ∈ Rn×n is defined as S ′ = T−1(S) := [ T−1AT T−1B ] . A linear transformation of a set V ⊆ Rn under the same mapping is defined as Y = T−1V := {y ∈ Rn | y = T−1v, v ∈ V}. 22 2.4. Approach I: Structure Decomposition Definition 2.5 (Unidirectionally Coupled). The LTI system that consists of two subsystems ẋ1 = A1x1 + ∆x2, (2.14) ẋ2 = A2x2 (2.15) with A1 ∈ Rk×k, A2 ∈ R(n−k)×(n−k), ∆ ∈ Rk×(n−k), x1(t) ∈ Rk, and x2(t) ∈ Rn−k, is said to be unidirectionally coupled since the trajectories of (2.14) are affected by those of (2.15), while (2.15) evolves independently from (2.14). The worst-case unidirectional coupling can thus be character- ized by the maximum row sum ‖∆‖. Definition 2.6 (Unidirectionally Weakly-Coupled). Let there be a non- singular transformation matrix T ∈ Rn×n, such that [zT1 , zT2 ]T = T−1[xT1 , xT2 ]T, and ż1 = A1z1 + ∆̃ z2 (2.16) ż2 = A2z2. (2.17) Then (2.16) and (2.17) are said to be unidirectionally weakly-coupled (in comparison to (2.14) and (2.15)) if ‖∆̃‖ ≤ ‖∆‖. (2.18) Definition 2.7 (Disjoint Input). Let there be a nonsingular transformation matrix T ∈ Rn×n and a coordinate space z = T−1x in which (2.12) can be partitioned into N subsystems as żi = Ãizi + gi(u), i = 1, . . . , N, (2.19) with gi(u) := B̃iu. The input u = [u1, . . . , up] T ∈ U ⊂ Rp is disjoint across these subsystems if ∀s ∈ {1, . . . , p}, ∀i, j ∈ {1, . . . , N}, i 6= j, ∂gi(u) ∂us 6= 0 → ∂gj(u) ∂us = 0. (2.20) 23 2.4. Approach I: Structure Decomposition Example 2.1. Consider the subsystem trajectories of żi = Ãizi + B̃iu with B̃i = [B̃i1, B̃i2, B̃i3], u = [u1, u2, u3] T, and i ∈ {1, 2}. The input u is disjoint if for example[ z1(t) z2(t) ] = [ eÃ1(t−t0) 0 0 eÃ2(t−t0) ][ z1(t0) z2(t0) ] + ∫ t t0 [ eÃ1(t−r) 0 0 eÃ2(t−r) ][ B̃11 0 0 0 B̃22 B̃23 ]u1(r)u2(r) u3(r) dr. (2.21) This ensures that no input from the input vector u occurs in both subsystems. On the other hand, the input in the following equation is non-disjoint since u2 influences the trajectories of both subsystems:[ z1(t) z2(t) ] = [ eÃ1(t−t0) 0 0 eÃ2(t−t0) ][ z1(t0) z2(t0) ] + ∫ t t0 [ eÃ1(t−r) 0 0 eÃ2(t−r) ][ B̃11 B̃12 0 0 B̃22 B̃23 ] u1(r)u2(r) u3(r) dr. (2.22) Definition 2.8 (ETUC). A subsystem is said to be externally-trivially- uncontrollable (ETUC) if it possesses a null input matrix. Finally, consider the following two lemmas. Lemma 2.1 ([123, Lem. 2.7]). The Sylvester equation EX +XF +H = 0, (2.23) with E ∈ Rk×k, F ∈ Rm×m, and H ∈ Rk×m, has a solution X ∈ Rk×m if and only if rank [ (FT ⊗ Ik) + (Im ⊗ E) − vec(H) ] = rank [ (FT ⊗ Ik) + (Im ⊗ E) ] , (2.24) 24 2.5. Approach II: Set-Theoretic Methods where ⊗ denotes the Kronecker product and vec(H) is a vector formed by stacking the columns of H below one another. This solution is unique if and only if the eigenvalue sum λi(E)+λj(F ) 6= 0, ∀i ∈ {1, ..., k}, ∀j ∈ {1, ...,m}. Lemma 2.2 (Real Schur Form [43, Thm’s 7.1.3 and 7.4.1], [114, 5R]). For any real matrix M ∈ Rn×n there exists an orthogonal matrix U ∈ Rn×n such that UTMU = M̃ is real upper quasi-triangular, and the eigenvalues of M are the eigenvalues of the block diagonals (each of dimension 2 or less) of M̃ . Furthermore, the matrix U can be chosen to order the eigenvalues arbitrarily. Remark 2.1. There always exists a partitioning of M̃ such that M̃ =[ M̃11 M̃12 0 M̃22 ] . The size of the partitions can be chosen as desired, so long as each block diagonal entry (maximum size 2× 2) of M̃ is completely covered by exactly one of the blocks on the diagonal of the partitioned matrix. 2.5 Approach II: Set-Theoretic Methods 2.5.1 Objective and Problem Statement For the generic, continuously-valued system (2.1) we define our second ob- jective as follows. Problem 2.2 (Objective II). Express the viability kernel V iab[0,τ ](K,U) in terms of the maximal reachable sets Reach]t(K,U), t ∈ [0, τ ] to enable the use of Lagrangian methods for the computation of the viability kernel and the synthesis of safety-preserving controllers. 2.5.2 Definitions and Preliminaries In this section we will present the definitions required for the treatment of the case in which (2.1) is a continuous-time system (cf. Section 5.1.1) of the form ẋ = f(x, u). (2.25) 25 2.6. Robust Reachability Analysis: Competing Inputs Definition 2.9. We say that a vector field f : X ×U → X is bounded on K if there exists a norm ‖·‖ : X → R+ and a real number M > 0 such that for all x ∈ K and u ∈ U we have ‖f(x, u)‖ ≤M . Definition 2.10. A partition P = {t0, t1, . . . , tn} of [0, τ ] is a set of distinct points t0, t1, . . . , tn ∈ [0, τ ] with t0 = 0, tn = τ and t0 < t1 < · · · < tn. Further, we denote • the number n of intervals [tk−1, tk] in P by |P |, • the size of the largest interval by ‖P‖ := max|P |k=1{tk+1 − tk}, and • the set of all partitions of [0, τ ] by P([0, τ ]). Definition 2.11. For a signal u : [0, τ ]→ U and a partition P = {t0, . . . , tn} of [0, τ ], define the tokenization of u corresponding to P as the set of func- tions {uk : [0, tk − tk−1]→ U}k such that uk(t) = u(t+ tk−1). (2.26) Conversely, for a set of functions {uk : [0, tk − tk−1]→ U}k, define their con- catenation u : [0, τ ]→ U as u(t) = uk(t− tk−1), t ∈ [tk−1, tk]. (2.27) Definition 2.12. The ‖·‖-distance of a point x ∈ X from a nonempty set A ⊂ X is defined as dist(x,A) := inf a∈A ‖x− a‖. (2.28) For a fixed set A, the map x 7→ dist(x,A) is continuous. 2.6 Robust Reachability Analysis: Systems with Competing Inputs While the majority of this thesis deals with systems with a single input, there are instances where we will make use of the differential game framework and 26 2.6. Robust Reachability Analysis: Competing Inputs assume adversarial inputs. For example, an “artificial” disturbance input is considered in Section 3.3 to perform a robust reachability analysis for one of the subsystems of the original deterministic system. Another example is when we extend the results presented in Chapter 5 to compute the robust version of the viability kernel (known as the discriminating kernel [18]) for systems with competing inputs in Chapter 6. Therefore in this section we will briefly lay out the definitions and preliminaries concerning systems with not only control input, but also uncertainties and/or external disturbances. Consider the continuous-time system ẋ(t) = f(x(t), u(t), v(t)), x(0) = x0 (2.29) with disturbance input v(t) ∈ V where V is a compact convex subset of Rmv . The vector field f : X × U × V → X is assumed to be Lipschitz in x and continuous in both u and v. Let U[0,t] be as in (2.2) and define V[0,t] := {v : [0, t]→ Rmv measurable, v(t) ∈ V a.e.}. (2.30) For every t ∈ [0, τ ], x0 ∈ X , u(·) ∈ U[0,t], and v(·) ∈ V[0,t], there ex- ists a unique trajectory xu,vx0 : [0, t] → X that satisfies the initial condition xu,vx0 (0) = x0 and the differential equation (2.29). As before, when clear from the context we shall drop the subscript and superscript from the trajectory notation. We assume that the disturbance input v is unknown but takes values on the (bounded) set V. Note that v ∈ V can also be used to capture any (unknown but bounded) uncertainties in the model. In a differential game setting the information pattern between the players (i.e. control input u and disturbance input v) is important and must be spec- ified. The control input follows a feedback strategy, i.e. u(t) = û(x(t), t) ∈ U ∀t ∈ [0, τ ]. We assume non-anticipative strategies for the disturbance in- put. A map ρ : U[0,t] → V[0,t] is non-anticipative for v if for every u(·), u′(·) ∈ U[0,t], u(s) = u ′(s) for a.e. s ∈ [0, t] implies ρ[u](s) = ρ[u′](s) for a.e. s ∈ [0, t] [27]. This results in a conservative formulation of our desired constructs by 27 2.6. Robust Reachability Analysis: Competing Inputs giving v a slight advantage over u. (cf. [89] for more detail.) Definition 2.13 (Robust Maximal Reachable Set). The robust maximal reachable set at time t is the set of initial states for which there exists a control for every disturbance such that the trajectories emanating from those states reach K exactly at time t: Reach]t(K,U ,V) := {x0 ∈ X | ∃u(·) ∈ U[0,t], ∀ρ : U[0,t] → V[0,t], xu,ρ[u]x0 (t) ∈ K}. (2.31) Definition 2.14 (Robust Minimal Reachable Tube). The robust minimal reachable tube over the horizon [0, τ ] is the set of initial states for which there exists a disturbance for every control such that the trajectories emanating from those states reach K at some time t ∈ [0, τ ]: Reach[[0,τ ](K,U ,V) := {x0 ∈ X | ∃ρ : U[0,τ ] → V[0,τ ], ∀u(·) ∈ U[0,τ ], ∃t ∈ [0, τ ], xu,ρ[u]x0 (t) ∈ K}. (2.32) Definition 2.15 (Discriminating Kernel). The (finite-horizon) discriminat- ing kernel5 of K is the set of all initial states in K for which there exists a control such that the trajectories emanating from those states remain within K for every disturbance for all time t ∈ [0, τ ]: Disc[0,τ ](K,U ,V) := { x0 ∈ X | ∃u(·) ∈ U[0,τ ], ∀ρ : U[0,τ ] → V[0,τ ], ∀t ∈ [0, τ ], xu,ρ[u]x0 (t) ∈ K } . (2.33) Note that with V = {0} the discriminating kernel reduces to the viability kernel under the deterministic system ẋ = f(x, u): Disc[0,τ ](K,U , {0}) ≡ V iab[0,τ ](K,U). (2.34) Finally, we will assume that the Isaac’s condition holds and therefore 5The infinite-horizon discriminating kernel DiscR+(K,U ,V) is also known as the robust maximal controlled-invariant subset [11]. 28 2.7. To Over- or Under-Approximate? the players’ order can be interchanged. As such, the discriminating kernel and the robust minimal reachable tube are duals of one another: Disc[0,τ ](Kc,U ,V) = (Reach[[0,τ ](K,U ,V))c. (2.35) 2.7 To Over- or Under-Approximate? Any approximations of the minimal reachable tube (or its robust version) must be an over-approximation since the target set K for this construct is generally deemed unsafe. The minimal reachable tube is the set of states that can become unsafe regardless of the control input applied. An under- approximation of this set would exclude states for which such property holds, falsely labeling them as safe. In contrast, the viability (or the discriminating) kernel must be under- approximated to ensure that every trajectory initiating from this set stays viable in K. The constraint set K for this construct is generally associated with safety. The states that belong to the viability kernel must retain the property that for each of them there exists at least one control policy that can keep the trajectory of the system in K. Over-approximating the viability kernel would include states for which this property does not hold, falsely labeling them as safe. The approximative techniques we present in this thesis will ensure that the minimal reachable tube is always over-approximated and that the via- bility kernel is guaranteed to be under-approximated. 29 Chapter 3 Schur-Based Structure Decomposition1 The first decomposition technique we present to address Problem 2.1 is in- spired by a model reduction algorithm for systems with unstable modes [82, 110]. Our method decomposes LTI systems into either completely de- coupled or weakly-coupled subsystems. Reachability analysis can be per- formed on each resulting subsystem independently. Back projecting and intersecting each of the lower-dimensional reachable tubes provides an over- approximation of the actual minimal reachable tube. A Sylvester equation (or an optimization problem) is solved in order to eliminate (or minimize) the coupling between the subsystems. Additional constraints are imposed when the control input is non-disjoint across subsystems, to prevent under- approximation of the (unsafe) minimal reachable tube. In addition, at the end of this section we will also provide conditions under which a subspace reachable tube remains unchanged for all time and show how this can be used in conjunction with the proposed Schur-based decomposition technique to yield an even further reduction of complexity for a class of systems. Outline Via Lemma 2.2, as in [104], we obtain an upper block triangu- lar A-matrix for (2.13). We then perform a second similarity transforma- tion and obtain a decoupled (or weakly-coupled) block diagonal matrix by solving a Sylvester equation (or an optimization problem). Therefore, we effectively decompose the system into two either completely decoupled or unidirectionally weakly-coupled subsystems. In the case where the decom- 1A version of this chapter has been published in [57, 58]. 30 3.1. Disjoint Control Input position is decoupled, the reachable tube is computed separately for each isolated subsystem. When the decomposed subsystems are unidirectionally weakly-coupled, the reachable tube is computed independently for the iso- lated subsystem, whereas for the remaining subsystem, the effect of coupling is accounted for by treating the coupling terms as disturbance and perform- ing reachability with competing inputs. For both decoupled and unidirec- tionally weakly-coupled decompositions, the intersection of back projections of the lower dimensional reachable tubes is an over-approximation of the ac- tual reachable tube in the transformed coordinate space. When the control input across the decomposed subsystems is non-disjoint, a constrained opti- mization problem is solved in order to make one of the subsystems ETUC. In the following analysis, we assume a partitioning of (2.13) that results in exactly two subsystems. However, the proposed method is generalizable to N subsystems by applying the same decomposition algorithm to each subsystem iteratively. A higher number of subsystems (i.e. iterated decom- position) may result in a more conservative over-approximation of the actual reachable tube. For k < n, we now apply Lemma 2.2 with transformation matrix U ∈ Rn×n to (2.13) to obtain [ Ã11 Ã12 B̃1 0 Ã22 B̃2 ] (3.1) with Ã11 ∈ Rk×k, Ã12 ∈ Rk×(n−k), Ã22 ∈ R(n−k)×(n−k), B̃1 ∈ Rk×p, and B̃2 ∈ R(n−k)×p. 3.1 Disjoint Control Input Consider the case in which the control input is disjoint across candidate subsystems. Proposition 3.1. If there exists a solution X ∈ Rk×(n−k) to the Sylvester equation Ã11X −XÃ22 + Ã12 = 0, (3.2) 31 3.1. Disjoint Control Input then a transformation W = [ Ik X 0 In−k ] ∈ Rn×n (3.3) makes (3.1) completely decoupled. Proof. cf. [82, 104, 110]. Applying the (invertible) transformation W to (3.1), we obtain[ Ã11 Ã11X−XÃ22+Ã12 B̂1 0 Ã22 B̂2 ] = [ Ã11 0 B̂1 0 Ã22 B̂2 ] . (3.4) Notice that the resulting subsystems [ Ã11 B̂1 ] and [ Ã22 B̂2 ] have been effectively decoupled through the coordinate transformation z = T−1x, T = UW . As we shall see in Section 3.3, reachability analysis (in this trans- formed coordinate space) can then be performed on each lower-dimensional subsystem separately. Now consider the case in which there is no solution to the Sylvester equation (3.2). Proposition 3.2. If (3.2) does not have a solution, then the transformation (3.3) with X = arg min Q∈Rk×(n−k) ‖Ã11Q−QÃ22 + Ã12‖ (3.5) results in unidirectionally weakly-coupled subsystems w.r.t. (3.1). Proof. Consider Ac := Ã11X −XÃ22 + Ã12 6= 0 in (3.4). It is clear that in the transformed coordinate space characterized by z = (UW )−1x, the state vector z2 evolves independently of z1 since ż2 = Ã22z2 + B̂2u2. However, z1 is affected by z2 through Ac. That is, we have ż1 = Ã11z1 + B̂1u1 + Acz2. (Note that ui, i ∈ {1, 2}, is the effective portion of the input vector u for the ith subsystem.) Minimization of the infinity norm of Ac therefore translates into minimizing, i.e. weakening, the worst-case unidirectional coupling of 32 3.2. Non-Disjoint Control Input z1 with z2. To see this, let X ∗ = arg min{‖Ã11Q − QÃ22 + Ã12‖ | Q ∈ Rk×(n−k)}. Then the hypothesis ‖Ã12‖ < ‖Ã11X∗ − X∗Ã22 + Ã12‖ would imply that X∗ = 0 can never be a solution. Since there are no constraints in (3.5) imposing this restriction, by contradiction we conclude that ‖Ã11X∗− X∗Ã22 + Ã12‖ ≤ ‖Ã12‖. Therefore, according to Definition 2.6, the resulting subsystems are unidirectionally weakly-coupled. Remark 3.1. The objective function of (3.5) is convex and therefore a solution always exists. Remark 3.2. The main rationale behind minimizing the infinity norm of the unidirectional coupling term (and thus, obtaining unidirectionally-weakly coupled subsystems) is that for the purpose of reachability analysis, the in- finity norm of this term will be used to formulate an upper-bound on the magnitude of the disturbance to the upper subsystem. This will be discussed further in Section 3.3. 3.2 Non-Disjoint Control Input Now consider a decomposition in which the control input is non-disjoint. In this case even if the dynamics of the subsystems are completely decoupled, their evolution is tightly paired through a common input. The difficulty arises, for example, when in the reachability computation a control value deemed optimal for one subsystem is in fact non-optimal for the full-order system. Blindly performing reachability for each subsystem separately may result in an under-approximation and additional measures have to be taken to ensure the over-approximation of the actual (unsafe) minimal reachable tube. One way to remedy this issue is by ensuring that at least one of the subsystems in the transformed coordinate space is ETUC. It is clear that in such a case the (otherwise non-disjoint) control action does not affect the evolution of the reachable tube of the ETUC subsystem. Therefore, an optimal control input for the subsystem with nonzero input matrix is also optimal for the full-order system. 33 3.3. Reachability in Lower Dimensions More formally, if either the pair (Ã22, B̂2) or the pair (Ã11, B̂1) in (3.4) is made ETUC, reachability analysis can be performed as in the disjoint control input case, separately for each subsystem. Assumption 3.1. C (B̃T1 ) ⊆ C (B̃T2 ) with C (·) the column-space operator. Proposition 3.3. The transformation (3.3) with X = arg min Q∈Rk×(n−k) ‖Ã11Q−QÃ22 + Ã12‖ (3.6) subject to QB̃2 = B̃1 results in unidirectionally coupled subsystems. Moreover, (Ã11, B̂1) is ETUC. Proof. Assumption 3.1 is the necessary and sufficient condition for solvabil- ity of the overdetermined equality constraint in (3.6) (cf. Appendix A.1). To see the trivial-uncontrollability of (Ã11, B̂1) consider B̂ := W −1B̃ in (3.4). We have [ B̂1 B̂2 ] := [ Ik −X 0 In−k ][ B̃1 B̃2 ] = [ B̃1 −XB̃2 B̃2 ] . (3.7) Constraining the optimizer in (3.6) to choose from the class of solutions{ X ∈ Rk×(n−k) | XB̃2 = B̃1 } simply enforces B̂1 = 0. The resulting subsystems can now be treated as in the disjoint control input case, and hence an over-approximation of the reachable tube in each subspace can be computed. 3.3 Reachability in Lower Dimensions Denote by S1 := Rk and S2 := Rn−k (3.8) the subspaces of Rn in which the two subsystems evolve. In the new coordi- nate space z = T−1x, T := UW reachability analysis can be performed on each lower-dimensional subsystem separately: 34 3.3. Reachability in Lower Dimensions Algorithm 3.1 Reachability in lower dimensions (Schur-Based) 1: Zτ ← T−1K 2: Z iτ ← ProjSi(Zτ ), ∀i ∈ {1, 2} . project onto ith subspace For lower subsystem: 3: Z2[0,τ ] ← Reach[[0,τ ](Z2τ ,U) For upper subsystem: 4: Treat Acz2 as disturbance . Ac := Ã11X −XÃ22 + Ã12 5: ζ ← ‖z2‖ ≡ supv∈Z2 [0,τ ] ‖v‖ 6: Compute upper-bound ‖Acz2‖ ≤ ‖Ac‖ζ 7: Z1[0,τ ] consrv.←− Reach[[0,τ ](Z1τ ,U ,B(‖Ac‖ζ)) 8: return(Z1[0,τ ],Z2[0,τ ]) Note that steps 4 through 6 of Algorithm 3.1 may or may not be needed depending on whether the subsystems are obtained from Propositions 3.1, 3.2, or 3.3. The following scenarios describe how the input(s) are quantified to construct the subsystem reachable tubes: S1 (Proposition 3.1 is used): For both Z1[0,τ ] and Z2[0,τ ], the single input is control and it is universally quantified. S2 (Proposition 3.2 is used): For Z1[0,τ ] the control input is universally quantified while the disturbance input (unidirectional coupling) is ex- istentially quantified. For Z2[0,τ ] the single input is control and it is universally quantified. S3 (Proposition 3.3 is used): For Z1[0,τ ] the single input is disturbance (unidirectional coupling) and it is existentially quantified. Note that in this case the robust minimal reachability operator in step 7 of the al- gorithm is effectively replaced by Reach][0,τ ](Z1τ ,B(‖Ac‖ζ)). For Z2[0,τ ] the single input is control and it is universally quantified. The over-approximation of the actual minimal reachable tube of the full- order system in X can be obtained using the following lemma. 35 3.3. Reachability in Lower Dimensions Lemma 3.1 ([91, 113]). Let Z i[0,τ ], i ∈ {1, 2}, be the computed lower- dimensional over-approximative reachable tube of subsystem i. Then the inverse transformation of the intersection of the back-projection of these sets onto Rn is a guaranteed over-approximation of the actual full-order reachable tube Reach[[0,τ ](K,U) of the system (2.13): T ( (Z1[0,τ ] × S2) ∩ (S1 ×Z2[0,τ ]) ) ⊇ Reach[[0,τ ](K,U). (3.9) 3.3.1 Formulating an Upper-Bound on Growth of Z1[0,τ ] in Scenario S3 When the subsystems are obtained via Proposition 3.3, the reachable tube in the subspace of the ETUC subsystem is computed without the need for solving a differential game. In fact, for this subsystem the unidirectional coupling is treated as disturbance and, therefore, it is existentially quanti- fied. Consequently, this disturbance together with the dynamics strive to enlarge the reachable (unsafe) set as much as possible. This allows us to for- mulate an analytic upper-bound on the over-approximation of the reachable tube in this subspace in terms of system and design parameters: Let z̄1 ∈ Z1τ and suppose D[0,τ ] is the set of measurable functions from [0, τ ] to B(‖Ac‖ζ). There exists an admissible input ϑ(·) ∈ D[0,τ ] such that (using time-reversed dynamics) we have z1 := e −Ã11τ z̄1 − ∫ τ 0 e−Ã11(τ−r)ϑ(r)dr, (3.10) z1 ∈ Z1[0,τ ]. (3.11) Bounding the effect of the input on the evolution of the trajectories we 36 3.3. Reachability in Lower Dimensions obtain ‖z1 − e−Ã11τ z̄1‖ ≤ ∫ τ 0 e‖Ã11‖(τ−r)‖Ac‖ζdr (3.12) = e‖Ã11‖τ − 1 ‖Ã11‖ ‖Ac‖ζ (3.13) = ( lim M→∞ M∑ i=1 τ i‖Ã11‖i−1 i! ) ‖Ac‖ζ (3.14) ≤ ( lim M→∞ M∑ i=1 τ i (√ k σ(Ã11) )i−1 i! ) ‖Ac‖ζ (3.15) =: µ[0,τ ] (3.16) where σ(·) is the largest singular value operator, and k is the dimension of the ETUC subsystem. Therefore, an upper-bound for how much Z1[0,τ ] can grow in backward time can be written as Z1[0,τ ] ⊆ ⋃ t∈[0,τ ] e−Ã11tZ1τ ⊕ B(µ[0,τ ]) (3.17) in which ⊕ denotes the Minkowski sum. In particular, the choice of k, the magnitude of the unidirectional coupling ‖Ac‖, the supremum of the reachable tube in the lower subspace ζ = supv∈Z2 [0,τ ] ‖v‖, and the largest singular value of the upper subsystem σ(Ã11) can all affect the conservatism of the reachable tube Z1[0,τ ]. Moreover, given k and τ , the flexibility of the Schur form in placing the eigenvalues in any order along the block-diagonals of Ã can be exploited to make this subsystem evolve with slower dynamics. Through various tests we were able to confirm that doing so could potentially prevent the excessive growth of Z1[0,τ ] by influencing both e−Ã11t and µ[0,τ ]. 37 3.4. Further Reduction of Complexity 3.4 Further Reduction of Complexity in Reachability for a Class of Unstable Systems We now demonstrate that for a specific class of unstable LTI systems, the Schur-based decomposition can be used to further reduce the computational burden associated with reachability analysis. Particularly, we decompose any full-order unstable system into stable and anti-stable subsystems with disjoint input across them. To do this, we employ the presented Schur-based decomposition while rearranging the order of eigenvalues such that the lower (controlled) subsystem contains only the non-negative eigenvalues and the upper (uncontrolled and possibly perturbed) subsystem contains the strictly-negative ones. As we will show in Proposition 3.4, under certain conditions, reachability analysis in the anti- stable subspace need not be performed since the target and the reachable tubes coincide for all time. Proposition 3.4. Suppose that for a controlled linear system (2.12) the following conditions are satisfied. (i) K is convex (but possibly arbitrarily shaped) and 0 ∈ K; (ii) the A-matrix is anti-stable (analytic in the open left-half complex plane) with repeated and real eigenvalues λ1 = · · · = λn ≥ 0; (iii) the algebraic and geometric multiplicities of λi(A) are equal. Then for any τ ∈ R+, Reach[[0,t](K,U) = K ∀t ∈ [0, τ ]. (3.18) Proof. The proof is provided in Appendix A.2. Remark 3.3. Condition (i) is easily generalizable to star-convex sets for which the origin is the convergence point (any line segment from the origin to x ∈ K is contained in K). An example of this is when the states are constrained to lp-space with 0 < p < 1. 38 3.4. Further Reduction of Complexity K nonconvex set λ λ= ji Figure 3.1: Phase-plane of a planar system with a non-convex target set K. Even though conditions (ii) and (iii) are satisfied, the backward reachable tube will grow. An intuitive 2-dimensional illustration of various cases that would vio- late conditions in Proposition 3.4 is given in Figures 3.1 and 3.2 where the trajectories are shown in backward time. Note that although Proposition 3.4 is stated in terms of a general full- order system (and as such, may seem too restrictive), it makes the following assertion: Corollary 3.1. If any isolated subsystem of any given unstable system in any coordinate space satisfies the conditions in Proposition 3.4, then the minimal reachable tube for that subsystem remains precisely equal to the target set in the respective subspace. Suppose that reachability analysis is to be performed for an unstable system ẋ = Ax + Bu, u ∈ U with k negative and (n − k) non-negative eigenvalues for a target set K. We apply Schur-based decomposition with an appropriately synthesized transformation matrix T to obtain[ Ã− Ac 0 0 Ã+ B̂2 ] (3.19) partitioned such that Ã+ and Ã− contain only non-negative and strictly- negative eigenvalues, respectively. If Ã+ and K satisfy the conditions (i), 39 3.5. Extension to Switched Linear Systems Re( ) 0 Im( ) 0 λ λ = ∧ =/ Im( ) 0 Re( ) 0 λ λ =/ ∧ =/ jiλ λ=/ jiλ λ= Algebraic Geometric>multiplicity multiplicity K Figure 3.2: Phase-plane with various eigenvalue scenarios that would violate conditions of Proposition 3.4 and thus causing the backward reachable tube to grow. (ii), and (iii), then according to Corollary 3.1 the reachable tube in the lower subspace does not grow and thus need not be computed. Reachabil- ity analysis is performed only for the upper subsystem resulting in further reduction of complexity by avoiding altogether the reachable tube compu- tation in the lower subspace. Specifically, step 3 in Algorithm 3.1 is entirely omitted. The over-approximation of the full-order reachable tube can then be calculated according to (3.9) with Z2[0,τ ] = Reach[[0,τ ](ProjS2(T−1K),U) ≡ ProjS2(T −1K) and Z1[0,τ ] = Reach][0,τ ](ProjS1(T−1K),B(‖Ac‖ζ)). Note that linear transformation preserves convexity. Therefore the pro- jection of the transformation of K onto the lower subspace (i.e. Z2τ := ProjS2(T −1K)) is convex if K is, and contains the origin if K does. 3.5 Extension to Switched Linear Systems The extension of our transformation-based method to switched dynami- cal systems is fairly straightforward. Consider the hybrid automaton H = (Q,X , f,U , E,R) with discrete modes Q = {qi}, continuous states x ∈ X , 40 3.6. Numerical Examples continuous control inputs u ∈ U , vector fields f : Q×X × U → X , (qi, x, u) 7→ Aix+Biu, (3.20) edges E ⊆ Q×Q, and (identity) reset maps R : E ×X → 2X . Let K(qi) and Reach[[0,τ ](qi,K(qi),U) denote the target and the reach- able tubes in mode qi ∈ Q, respectively. Also let Ti be the transformation matrix for mode qi obtained from the Schur-based decomposition technique described previously. For simplicity of presentation we assume that H has only two modes qi and qj such that (qi, qj) ∈ E. The backward reachable tube in mode qj can be directly expressed as TjReach [ [0,τ ] ( qj , T −1 j TiReach [ [0,τ ] ( qi, T −1 i K(qi),U ) ,U ) . (3.21) Reachability analysis is then performed on lower-dimensional subsystems in each mode according to Algorithm 3.1. 3.6 Numerical Examples Although complexity reduction through Schur-based decomposition can be used in conjunction with any reachability/viability technique that can ac- commodate both existentially and universally quantified inputs, we demon- strate the applicability and practicality of our method using a number of examples (up to 8D) that employ the Level Set Toolbox (LS) [90]. While LS has mainly been used for systems of low dimensionality [8], our complex- ity reduction approach can facilitate the use of LS for higher dimensional systems for which safety-preserving controller synthesis and/or handling of non-convex, arbitrarily-shaped sets is important. All computations are performed on a dual core Intel-based computer with 2.8 GHz CPU, 6 MB of cache and 3 GB of RAM running single-threaded 32-bit Matlab 7.5. 41 3.6. Numerical Examples Figure 3.3: The non-convex target set Zτ in the transformed coordinate space. 3.6.1 Arbitrary 3D System Consider an arbitrary 3D LTI system with A = −0.5672 −0.7588 −0.62823.1364 −1.1705 2.3247 1.8134 −1.7689 −2.6930 , B = 0.0731 −0.1639−0.7377 −0.3578 0.1470 0.2410 and input u = [u1, u2] T ∈ R2, ‖u‖ ≤ 1.1. We choose a non-convex target (unsafe) set K ⊂ R3 such that in the transformed coordinate space we have Zτ = T−1K as shown in Figure 3.3. Here, T is the transformation matrix obtained through Proposition 3.1 that decomposes the system into two subsystems (one 2D and one 1D) with disjoint control across them: T−1AT = −1.6653 −3.4560 0 1.8706 −1.4653 0 0 0 −1.3000 , T−1B = −0.7530 0 0.0640 0 0 0.2500 . Hence, the decoupled subsystems are ż1 = [−1.6653 −3.4560 1.8706 −1.4653 ] z1 + [−0.7530 0.0640 ] u1 and ż2 = [−1.3000 ] z2 + [ 0.2500 ]u2. We obtain an over-approximation of the actual reachable tube, as shown in Figure 3.4. Reachability calculation is performed over a grid with 101 42 3.6. Numerical Examples (a) Comparison in z-space (b) (z1, z2) cross-section (c) (z1, z3) cross-section (d) (z2, z3) cross-section Figure 3.4: Schur-based over-approximation (transparent light) vs. actual (solid dark) reachable tubes in the transformed coordinate space for Example 3.6.1. nodes in each dimension for τ = 2 s. The computation time for the actual and the Schur-based reachable tubes (including decomposition and projec- tions) were 5823.73 s and 22.87 s, respectively. 43 3.6. Numerical Examples 3.6.2 4D Aircraft Dynamics Consider longitudinal aircraft dynamics ẋ = Ax+Bδe, A = −0.0030 0.0390 0 −0.3220 −0.0650 −0.3190 7.7400 0 0.0200 −0.1010 −0.4290 0 0 0 1 0 , B = 0.0100 −0.1800 −1.1600 0 with state x = [u, v, θ̇, θ]T ∈ R4 comprised of deviations in aircraft velocity [ft/s] along and perpendicular to body axis, pitch-rate [crad/s], and pitch angle [crad] respectively2, and with input δe ∈ [−13.3◦, 13.3◦] ⊆ R the ele- vator deflection. These matrices represent stability derivatives of a Boeing 747 cruising at an altitude of 40 kft with speed 774 ft/s [16]. Define a non-convex target (unsafe) set K such that in the transformed coordinate space Zτ = {z ∈ R4 | ‖z‖ > 0.15, z = T−1x, x ∈ K} where T is the transformation matrix obtained through our method. We first de- compose the system into two 2D subsystems. Since the control input is non-disjoint across the resulting subsystems, we use Proposition 3.3 and obtain unidirectionally coupled subsystems, one of which is ETUC (see Ap- pendix A.3). The reachability calculation is performed over a grid with 41 nodes in each dimension for τ = 5 s. The computation time for the actual and the Schur-based minimal reachable tubes (including decomposition and projections) were 28546.80 s and 54.64 s, respectively—a significant reduc- tion. Since the computed sets are 4D, we plot a series of 3D snapshots of these 4D objects at specific values of z4 (Figure 3.5). The aircraft flight envelope (safe) is represented by the area inside of the shaded regions. 3.6.3 8D Distillation Column Consider the dynamic model of a binary distillation column obtained from [111] with 2crad = 0.01 rad ≈ 0.57◦. 44 3.6. Numerical Examples Figure 3.5: Schur-based (solid dark) vs. actual (transparent light) finite- horizon viability kernels (safe) in the transformed coordinate space for Example 3.6.2. The computed minimal reachable tube and its over-approximation are the non-convex complements of these objects. A = −0.5774 3.0567 0.0073 −0.8121 0.3034 −0.3035 0.0072 −0.1542 −2.7290 −0.7147 −0.3430 1.5321 0.6643 0.2896 −0.0013 0.0926 0 0 −0.3891 −0.9956 0.0182 0.0235 0.0049 0.0506 0 0 1.3640 −1.3363 −0.9037 −0.4686 −0.0009 −0.1887 0 0 0 0 −0.7357 −0.2275 −0.0082 −0.0021 0 0 0 0 0 −0.2259 0.0021 −0.0457 0 0 0 0 0 0 −0.0052 0.0024 0 0 0 0 0 0 0 −0.0755 B = [ −0.0335 −0.4534 −0.8005 0.5497 1.2886 0.3132 0.7117 0.0599 −0.1228 −0.0711 −0.2612 −0.1344 −0.0504 −0.2249 −0.6994 −0.3014 ]T . The input u = [u1, u2] T ∈ R2 with u1, u2 ∈ [0, 1] is comprised of reflux and boilup flows, respectively. The full-order system with state vector x ∈ 45 3.6. Numerical Examples R8 is first decomposed into two (unidirectionally coupled) 4D subsystems using Proposition 3.3, since the control vector is non-disjoint across the two candidate subsystems. Similarly, each of these 4D subsystems is decomposed into two 2D subsystems. Since the upper 4D subsystem is made ETUC through (3.6), its decomposition is disjoint and therefore Proposition 3.1 is used to obtain the 1st and 2nd (decoupled) 2D subsystems. On the other hand, for the lower 4D subsystem the decomposition results in non-disjoint control input. Therefore Proposition 3.3 is employed and the 3rd and 4th (unidirectionally coupled) 2D subsystems are obtained (see Appendix A.3). Reachability is first performed on the 3rd and 4th subsystems while taking the effect of unidirectional coupling into account. Next, the reachable tubes of the 1st and 2nd subsystems are computed while treating the effect of the 3rd and 4th subsystems as disturbance. We label the 2D transformed state subspaces as w̃1 = [w1, w2] T, w̃2 = [w3, w4] T, q̃1 = [q1, q2] T, and q̃2 = [q3, q4] T. Notice that R4 3 q = [q̃T1 , q̃T2 ]T = T−13 z̃2, R4 3 w = [w̃T1 , w̃T2 ]T = T−12 z̃1, and R8 3 z = [z̃T1 , z̃T2 ]T = T−11 x with z̃1, z̃2 ∈ R4. We assume that the target (unsafe) set K ⊂ R8 is chosen such that the transformations T−11 ∈ R8×8, T−12 ∈ R4×4, and T−13 ∈ R4×4 result in Wτ := {w ∈ R4 | ‖w‖ > 20} and Qτ := {q ∈ R4 | ‖q‖ > 20}. The target sets for the 2D subsystems is simply the projection ofWτ and Qτ onto their corresponding subspaces. Lower dimensional reachability is performed over a grid with 101 nodes in each dimension for τ = 6 s. The overall computation time (including decomposition and projections) was 94.31 s. The complement of the shaded regions in Figure 3.6 over-approximate the reachable (unsafe) set in each of the 2D subspaces. The full 8D reachable tube is the intersection of the back-projection of the 2D reachable tubes. The actual (minimal) reachable tube is not shown since it is prohibitively computationally expensive to compute using any Eulerian method including LS. 46 3.6. Numerical Examples q3 q 4 −20 0 20 −20 0 20 q1 q 2 −20 0 20 −20 0 20 w3 w 4 −20 0 20 −20 0 20 w1 w 2 −20 0 20 −20 0 20 Figure 3.6: The Schur-based viability kernel (safe) of Example 3.6.3 in trans- formed 2D subspaces. 3.6.4 4D Unstable System (An Example for Section 3.4) Consider an unstable system [51, Ex. 2.2.1] with A = 0 1 0 0 0.1023 0 −0.0085 0 0 0 0 1 −0.0153 ε 0.0993 0 , B = 10−3× 0 −0.8696 0 0.1304 . Let the eigenvalues of the system be slightly perturbed as determined by parameter ε ∈ R. With ε = 0.0491 the real anti-stable eigenvalues coincide. We apply the method described in Section 3.4 and obtain two 2D sub- systems (with separated stable and anti-stable eigenvalues) across which the 47 3.7. Summary and Further Discussions input is disjoint. The system matrices in the transformed coordinates are T−1AT= −0.3426 0.0354 −0.6988 0.1399 −0.0000 −0.2912 0.9481 −0.0000 0 0 0.3150 −0.0135 0 0 0.0003 0.3188 , T−1B=10−3× 0 0 −0.7621 0.3426 . A target (unsafe) set K is chosen such that Zτ = {z ∈ R4 | √ zTz ≤ 0.2, z = T−1x, x ∈ K}, i.e. a small Euclidean ball of radius 0.2 around the origin. The magnitude of the input is bounded by |u| ≤ 1. Using reachability analysis we attempt to identify the set of initial states that reach Zτ in τ = 3 seconds, regardless of the input applied. Since all conditions in Proposition 3.4 are satisfied for the lower sub- system, to obtain an over-approximation of the full-order system, we only compute the over-approximation of the reachable tube in its stable sub- space. The minimal reachable tube and its over-approximation are shown in Figure 3.7. Reachability was performed over a grid with 41 nodes in each dimension. The overall computation time (including decomposition and projections) was 2.8 s. In comparison, computing the reachable tube of the full-order system would require 1741.6 s. 3.7 Summary and Further Discussions In this chapter we presented our first decomposition technique, Schur-based decomposition, to facilitate a comparatively more scalable computation of the minimal reachable tube (and by duality the viability kernel) for LTI systems using Eulerian methods. The decomposition was evaluated in terms of whether the resulting sub- systems had disjoint or non-disjoint control inputs. In the event that a Sylvester equation can be solved, the decomposition yields two decoupled subsystems. When the Sylvester equation cannot be solved, its infinity norm minimization yields two weakly coupled subsystems. Additional constraints are considered for the case in which the control input is non-disjoint across 48 3.7. Summary and Further Discussions Figure 3.7: 3D snapshots of the actual (solid dark) unsafe reachable tube vs. its over-approximation (transparent light) in the transformed coordinate space for Example 3.6.4. The over-approximation was computed using Schur-based decomposition in conjunction with Proposition 3.4 for only one of the subsystems. decomposed subsystems. Reachability analysis is then performed on the re- sulting subsystems independently. We applied this technique to a variety of examples computed with the Level-Set Toolbox, and found computational time significantly reduced when our method was employed. Furthermore, we presented conditions under which the minimal reachable tube and the target set coincide. We then showed that the proposed Schur-based decom- position can be used together with these conditions in order to significantly reduce the computational complexity of reachability analysis for a class of unstable systems. The introduced conservatism in the over-approximation of the minimal reachable tube can be mitigated to some degree by considering a time- dependent formulation of the disturbance to the upper subsystem and per- forming reachability in sub-intervals of [0, τ ]. This procedure is described 49 3.7. Summary and Further Discussions towards the end of the next chapter for our second decomposition technique. 50 Chapter 4 Riccati-Based Structure Decomposition1 Our second proposed decomposition technique that aims to address Prob- lem 2.1 draws upon the so-called Riccati transformation—a two-stage coor- dinate transformation based on the solutions of a nonsymmetric algebraic Riccati equation (NARE) and a Sylvester equation. This transformation, originally introduced in [19] for decoupling of singularly perturbed systems, was later generalized in [62] to larger classes of autonomous LTI systems. An in-depth overview of the application of this transformation in optimal control theory, singular perturbation theory, and asymptotic approximation theory can be found in [112], while more recent advances are given in [32] and [107]. Outline When the transformation results in input that is disjoint across the candidate subsystems, the standard Riccati transformation can be used to decompose the system into two decoupled subsystems. For the case in which the control input is non-disjoint across the decomposed subsystems, we propose a modified Riccati transformation (an extension to the stan- dard Riccati transformation) which, in addition to parameterizing the uni- directional coupling between the subsystems, makes one of the subsystems ETUC. In the new coordinate space reachability analysis can then be per- formed in lower dimensions for each subsystem separately. The intersection of back projections of the lower dimensional reachable tubes is an over- approximation of the actual reachable tube in the transformed coordinate 1A version of this chapter has been published in [56]. 51 Chapter 4. Riccati-Based Structure Decomposition space. In the following analysis we assume a partitioning of (2.13) that results in exactly two subsystems. However, the proposed method can be generalized to N subsystems by applying the same decomposition algorithm iteratively (see Section 4.3). Let (2.13) be partitioned as S = [ A11 A12 B1 A21 A22 B2 ] (4.1) with A11 ∈ Rk×k, A12 ∈ Rk×(n−k), A21 ∈ R(n−k)×k, A22 ∈ R(n−k)×(n−k), B1 ∈ Rk×p, and B2 ∈ R(n−k)×p, for some k < n. Now consider the nonsin- gular transformation matrices T1 = [ Ik 0 −L In−k ] ∈ Rn×n, (4.2) T2 = [ Ik M 0 In−k ] ∈ Rn×n. (4.3) With L ∈ R(n−k)×k and M ∈ Rk×(n−k) that satisfy (NARE:) R(L) := LA11 −A22L− LA12L+A21 = 0, (4.4) (Sylvester:) S (M) := ( A11 −A12L ) M −M(A22 + LA12)+A12 = 0, (4.5) the transformed system is S ′ = T−11 (S) = A11 −A12L A12 B1 *0 R(L) A22 + LA12 LB1 +B2 , (4.6) S ′′ = T−12 (S ′) = [ A11 −A12L :0S (M) (I −ML)B1 −MB2 0 A22 + LA12 LB1 +B2 ] . (4.7) 52 4.1. Disjoint Control Input Solutions to (4.4) and (4.5) may not always exist. 4.1 Disjoint Control Input Consider the resulting subsystems S ′′1 = [ A11 −A12L (I −ML)B1 −MB2 ] , (4.8) S ′′2 = [ A22 + LA12 LB1 +B2 ] . (4.9) Suppose the control input is disjoint across these dynamically decoupled can- didate subsystems (that is, no common input occurs in both subsystems). If the following assumption regarding the input set is satisfied, reachability analysis for each subsystem can be performed independently. Paralleliza- tion of reachability calculations in each subspace could further reduce the computational time. Assumption 4.1. U = ∏2i=1 Ui where Ui is the subset of Rp from which the portion of the input vector u acting on subsystem i draws its values. Assumption 4.1 enures that the inputs acting on the two subsystems are independent of one another. Note that this condition is satisfied for most physical systems where actuators are commonly uncorrelated. In the general case, however, we can under-approximate U by a hyper-rectangle formed by the direct product of p one-dimensional intervals. The resulting reachable tube in each subspace will be a conservative over-approximation of the actual reachable tube in that subspace since for any given system, a smaller input authority yields a larger (minimal) reachable tube. 4.2 Non-Disjoint Control Input When the control input across the candidate subsystems is non-disjoint, even if the states of the subsystems are completely decoupled, their evo- lution is tightly paired through common input. Difficulty arises, for ex- ample, when in the reachability computation a control value deemed op- 53 4.2. Non-Disjoint Control Input timal for one subsystem is in fact non-optimal for the full-order system. Blindly performing reachability for each subsystem separately may result in an under-approximation, so appropriate measures must be taken to ensure over-approximation of the actual (unsafe) reachable tube. One way to remedy this issue is by ensuring that at least one of the subsystems in the transformed coordinate space is ETUC. It is clear that in such a case the (otherwise non-disjoint) control action does not affect the evolution of the ETUC subsystem. Therefore, an optimal input for the subsystem with nonzero input matrix is also optimal for the full-order sys- tem. We propose a modified Riccati transformation that ensures that one of the subsystems in the transformed coordinates is ETUC, hence reachability analysis can be performed separately for each subsystem in its corresponding lower dimensional subspace. Remark 4.1. With an ETUC subsystem, Assumption 4.1 on the shape of the input set is lifted: The input u acts only on one of the subsystems, therefore the shape of U becomes irrelevant. 4.2.1 Transformation 1 (ETUC Subsystem) Consider a transformation through which the lower subsystem can be made ETUC. That is, in (4.6) for the transformation matrix T1 we seek an L in R(L) that is also a solution of LB1 +B2 = 0. Assumption 4.2. C (BT2 ) ⊆ C (BT1 ) with C (·) the column-space operator. Lemma 4.1. Under Assumption 4.2, the class of solutions of LB1 = −B2 w.r.t. L ∈ R(n−k)×k can be characterized by L := { −B2B†1 + Z − ZB1B†1, Z ∈ R(n−k)×k } , (4.10) where † denotes the Moore-Penrose pseudoinverse. Proof. cf. [103] and [45]. Assumption 4.2 is the necessary and sufficient condition for solvability of LB1 = −B2. (See Appendix A.1.) 54 4.2. Non-Disjoint Control Input Substituting (4.10) for L in R(L) we obtain R̂(Z) := ZΠ + Γ + Z ( A12 −B1B†1A12 ) Z(B1B † 1 − I) + ( A22 −B2B†1A12 ) Z(B1B † 1 − I), (4.11) where Π = −(B1B†1 − I) ( A11 +A12B2B † 1 ) , (4.12) Γ = ( A22B2B † 1 +A21 ) −B2B†1 ( A12B2B † 1 +A11 ) . (4.13) To eliminate the non-invertible term (B1B † 1−I) from the right-hand side of (4.11) we equate R̂(Z) to some rank correcting term δF (Z) with F (Z) := Z ( A12 −B1B†1A12 ) Z + ( A22 −B2B†1A12 ) Z (4.14) and δ ∈ R\{−1, 0} a finite (but possibly large) parameter such that (B1B†1− (δ + 1)I ) is nonsingular: R̂(Z) = ZΠ + Γ + Z ( A12 −B1B†1A12 ) Z(B1B † 1 − I) + ( A22 −B2B†1A12 ) Z(B1B † 1 − I) (4.15) = ZΠ + Γ +F (Z)(B1B † 1 − I) .= δF (Z). (4.16) Simple algebraic manipulation and post-multiplication of R̂(Z)−δF (Z) = 0 by ( B1B † 1 − (δ + 1)I )−1 then results in a NARE in the variable Z: R1(Z) := ZÃ11 − Ã22Z − ZÃ12Z + Ã21 = 0 (4.17) with Ã11 = Π ( B1B † 1 − (δ + 1)I )−1 , Ã21 = Γ ( B1B † 1 − (δ + 1)I )−1 , Ã12 =( B1B † 1A12 −A12 ) , and Ã22 = ( B2B † 1A12 −A22 ) . Proposition 4.1. If a root Z ∈ R(n−k)×k of the NARE (4.17) exists, it 55 4.2. Non-Disjoint Control Input constitutes an L ∈ L that simultaneously satisfies LB1 +B2 = 0, (4.18a) R(L) = LA11 −A22L− LA12L+A21 = δF (Z). (4.18b) Proof. By virtue of (4.16), a matrix Z that satisfies (4.17) also satisfies (4.18) via (4.10). Remark 4.2. If p ≥ k in partitioning of (4.1), the set L reduces to the singleton {−B2B†1} and the method still applies. Theorem 4.1. The transformation (4.2) with L ∈ R(n−k)×k obtained through Proposition 4.1 makes the lower subsystem in (4.1) ETUC. Moreover, the coupling terms are altered such that the effect of the upper subsystem on the evolution of the lower subsystem is parameterized by δ. Proof. S ′ = T−11 (S) = [ A11 −A12L A12 B1 LA11 −A22L− LA12L+A21 A22 + LA12 LB1 +B2 ] (4.19) = [ A11 −A12L A12 B1 δF (Z) A22 + LA12 0 ] . (4.20) Remark 4.3. Note that the imposed δ-parameterization of the off-diagonal term δF (Z) in (4.20) provides an additional degree of freedom in adjusting (minimizing) the coupling of the two subsystems in the new coordinates. This will be discussed further in Section 4.2.3. Nonsymmetric Riccati equations have long been an active area of re- search [31]. To solve (4.17) we draw on the fixed-point algorithm described in [62] and derive the necessary conditions for the existence and uniqueness of a real root Z. 56 4.2. Non-Disjoint Control Input Suppose ( B2B † 1A12 −A22 ) is invertible. Define initial values as Z0 := ( B2B † 1A12 −A22 )−1 Γ ( B1B † 1 − (δ + 1)I )−1 , (4.21) A0 := Π ( B1B † 1 − (δ + 1)I )−1 − (B1B†1A12 −A12)Z0. (4.22) To find Z we look for D := Z − Z0 (4.23) by solving R̃1(D) := DA0 − ( B2B † 1A12−A22 + Z0 ( B1B † 1A12 −A12 )) D −D(B1B†1A12 −A12)D + Z0A0 = 0. (4.24) Lemma 4.2 ([62, Lem. 1]). Suppose ( B2B † 1A12 −A22 ) is nonsingular. If ∥∥(B2B†1A12 −A22)−1∥∥ ≤ 1 3 ( ‖A0‖+ ‖B1B†1A12 −A12‖‖Z0‖ ) (4.25) then (4.24) has a unique real root D that satisfies 0 ≤ ‖D‖ ≤ 2‖A0‖‖Z0‖ ‖A0‖+ ‖B1B†1A12 −A12‖‖Z0‖ (4.26) and is the fixed-point solution of the contraction Dk+1 = P1(Dk) given by P1(Dk) := ( B2B † 1A12 −A22 )−1( Z0A0 +DkA0 − Z0 ( B1B † 1A12 −A12 ) Dk −Dk ( B1B † 1A12 −A12 ) Dk ) . (4.27) Remark 4.4. Similarly to [62], it can be shown that the relative error ek := ‖Dk −D‖/‖D‖ after k iterations is bounded above by ek ≤ ( 3 ∥∥(B2B†1A12 −A22)−1∥∥(‖A0‖+ ‖B1B†1A12 −A12‖‖Z0‖))k (4.28) 57 4.2. Non-Disjoint Control Input and decreases as |δ| increases since ‖A0‖ and ‖Z0‖ are inversely related to |δ|. For a given δ, using D0 = 0 as initial condition we compute D iteratively. The fixed-point solution D∗ = P1(D∗) is then used to obtain Z = D∗ + Z0 which in turn solves R1(Z) = 0 in (4.17) and results in a matrix L, through (4.10), that satisfies both equations in (4.18). 4.2.2 Transformation 2 (Unidirectionally Coupled Subsystems) Consider the NARE R2(M) = ( A11 −A12L ) M −M(A22 + LA12) −M(δF (Z))M +A12 = 0. (4.29) For a given L, δ, and Z, if there exists a solution M that satisfies (4.29), we obtain the following: Theorem 4.2. The transformation (4.3) with M ∈ Rk×(n−k) satisfying NARE (4.29) makes the subsystems in (4.20) unidirectionally coupled. Proof. S ′′ = T−12 (S ′) = [ A11 −A12L−MδF (Z) :0R2(M) B1 δF (Z) A22 + LA12 + δF (Z)M 0 ] . (4.30) Remark 4.5. In the transformed coordinates the lower subsystem remains ETUC. Furthermore, the δ-parameterization of the unidirectional coupling between subsystems is also preserved. Before further analyzing the unidirectional coupling term δF (Z), let us derive the necessary conditions for the existence and uniqueness of a 58 4.2. Non-Disjoint Control Input solution M to (4.29) to be used with the same convergent iterative procedure described previously. For a given δ, Z, and L, let ( A11 − A12L ) be invertible and the initial values be defined as M0 := − ( A11 −A12L )−1 A12, (4.31) N0 := A22 + LA12 + δF (Z)M0. (4.32) We seek M by forming J := M −M0 (4.33) and solving R̃2(J) := JN0− ( A11−A12L−δM0F (Z) ) J+δJF (Z)J+M0N0 = 0. (4.34) Lemma 4.3 ([62, Lem. 1]). Suppose ( A11 −A12L ) is nonsingular. If ∥∥(A11 −A12L)−1∥∥ ≤ 1 3 ( ‖N0‖+ ‖δF (Z)‖‖M0‖ ) (4.35) then (4.34) has a unique real root J that satisfies 0 ≤ ‖J‖ ≤ 2‖N0‖‖M0‖‖N0‖+ ‖δF (Z)‖‖M0‖ (4.36) and is the fixed-point solution of the contraction Jk+1 = P2(Jk) given by P2(Jk) := ( A11 −A12L )−1( M0N0 + JkN0 + δM0F (Z)Jk + δJkF (Z)Jk ) . (4.37) Remark 4.6. As in [62], we can show that the relative error ek := ‖Jk − J‖/‖J‖ after k iterations is bounded above by ek ≤ ( 3 ∥∥(A11 −A12L)−1∥∥(‖N0‖+ ‖δF (Z)‖‖M0‖))k (4.38) and decreases as ‖δF (Z)‖, ‖A22‖, and ∥∥(A11 − A12L)−1∥∥ decrease. This 59 4.2. Non-Disjoint Control Input occurs when the ill-conditioning of the A-matrix increases (e.g. in the case of two-time-scale systems; see [63] and the references therein) and δ is chosen such that ‖δF (Z)‖ is minimized. Using J0 = 0 as initial condition we compute J iteratively. The fixed- point solution J∗ = P2(J∗) is then used to obtain M = J∗ + M0 which in turn solves R2(M) = 0 in (4.29). Note that both conditions (4.25) and (4.35) are conservative and their satisfaction ensures rapid convergence (usually within 2 or 3 iterations). In practice, the right-hand-side of these inequalities can be relaxed up to 10 times in most cases without causing divergence. 4.2.3 The Unidirectional Coupling Term (Choosing δ) Finally, we analyze the unidirectional coupling term δF (Z) and its behavior with respect to the free parameter δ. Since Z is an implicit function of δ, we adopt the extended notation δF (Z(δ)) to reflect this dependency. First, we formalize a conservative upper-bound on ‖δF (Z(δ))‖ as an explicit function of δ. This assures that the unidirectional coupling remains bounded for almost all admissible values of the free parameter δ. Proposition 4.2. The worst-case unidirectional coupling between the two subsystems in the transformed coordinates, i.e. ‖δF (Z(δ))‖ in (4.30), is (conservatively) bounded above such that ‖δF (Z(δ))‖ ≤ 1|δ| ( |δ|+ 1 |δ + 1| )2 a+ ( |δ|+ 1 |δ + 1| ) b ∀δ ∈ R\{−1, 0}, (4.39) where the constants a and b are independent of δ and are determined by a := α(b/β)2, b := 3‖B1B†1‖γβ, γ := ‖Γ‖ ∥∥(A22 − B2B†1A12)−1∥∥, α := ‖A12 −B1B†1A12‖, and β := ‖A22 −B2B†1A12‖. Proof. The proof is provided in Appendix B.1. Now consider inequalities (4.25) and (4.35), which are dependent on δ. Adequately chosen and sufficiently large values of δ help ensure that these 60 4.2. Non-Disjoint Control Input conditions are met. On the other hand, choosing δ exceedingly large de- feats the purpose of δ-parameterization of the unidirectional coupling term, since it can be shown that as δ grows, ‖δF (Z(δ))‖ approaches a problem- dependent constant that may not necessarily be an extremum point. Proposition 4.3. lim δ→±∞ ‖δF (Z(δ))‖ = ‖Γ‖ with Γ given by (4.13). Proof. This proof is also provided in Appendix B.1. It follows from Proposition 4.3 that 0 ≤ infδ‖δF (Z(δ))‖ ≤ ‖Γ‖. There- fore naively letting |δ| → ∞ essentially removes the added flexibility as- sociated with the δ-parameterization in the modified Riccati approach and instead enforces a trivial solution L = −B2B†1. While for some systems this solution may yield the smallest possible unidirectional coupling between the resulting subsystems (i.e. a unidirectional coupling with the lowest infinity norm), in most cases a carefully chosen δ not only facilitates the satisfaction of the convergence conditions (4.25) and (4.35), but also further minimizes the worst-case unidirectional coupling. Thus, formulated as an optimization problem, we seek a δ that solves the following: minimize δ∈R\{−1,0} f(δ) := ‖δF (Z(δ))‖ subject to (4.25) and (4.35). Note that this is a non-convex problem, and in general, f(·) may be a non-smooth function of δ. However, a global optimum need not be com- puted. Any suboptimal solution can be used as long as that solution yields a satisfactory degree of unidirectional coupling between the subsystems in the transformed coordinates. In addition, an approximation to the optimum point can be obtained numerically, for example by fine-griding the real line or using the bisection algorithm. In practice, while the exact shape of the function f(·) is problem-dependent, we have found (but have not proven) that in most cases it exhibits a be- havior similar to that of an absolute value proper rational function (over a 61 4.3. Recursive Decomposition −1e6 −1e4 −1e2 1e2 1e4 1e6 2 2.5 3 3.5 4 in fe a si b le re g io n δ f (δ ), f̂ (δ ) ||δF (Z(δ))|| ||Γ|| approximation Figure 4.1: The worst-case unidirectional coupling f(δ) = ‖δF (Z(δ))‖ (×’s) and its approximation f̂(δ) = |−27.65δ + 0.55| + 1.82 (dashed) computed for Example 4.5.1. The interval (−15,+15) over which (4.25) and (4.35) are violated is labeled as “infeasible region”. The asymptote limδ→±∞ f(δ) = ‖Γ‖ (dash-dotted) is also shown. The minimum of f(δ) occurs when δ ≈ +50. discontinuous domain) of the form f̂(δ) = ∣∣∣ c0 δk + c1 ∣∣∣+ c2 ∀δ ∈ D, (4.40) where D ⊂ R\ {−1, 0} is the union of the two segments of the real line for which the magnitude of δ is large enough such that (4.25) and (4.35) are both satisfied, k ∈ N, k : odd, c0 = −c1(δ∗)k, δ∗ = arg minδ∈Y f(δ), c2 = minδ∈Y f(δ), and c1 = ( limδ→±∞ f(δ) )− c2 = ‖Γ‖ − c2. For example, consider the transformed system in Section 4.5.1. Figure 4.1 shows f(δ) and its approximation f̂(δ) = |−27.65δ + 0.55|+ 1.82 evaluated where (4.25) and (4.35) hold. 4.3 Recursive Decomposition To apply the described decomposition technique in a recursive fashion, con- sider the resulting subsystems in (4.8) and (4.30). A recursive decomposition when the standard Riccati transformation can be used is straightforward. 62 4.4. Reachability in Lower Dimensions Suppose that the modified Riccati transformation is used throughout the process. In deeper level recursions, the decomposition can be applied to the uppermost subsystem since that subsystem is controlled whereas every other subsystem is ETUC. For example, to decompose a 6D system into three 2D subsystems, in the first recursion level, the partitioning can be chosen such that the resulting upper (controlled) subsystem is 4D and the lower (ETUC) subsystem is 2D. In the second recursion level, if the solutions exist, the 4D subsystem is then decomposed into two 2D subsystems. Note that in the recursive application of the decomposition, when the modified Riccati transformation is employed, all subsystems but one are ETUC. Therefore, this iterated decomposition may result in a more conser- vative over-approximation of the actual reachable tube. 4.4 Reachability in Lower Dimensions We will focus mainly on over-approximating the minimal reachable tube that can only be computed using the computationally intensive Eulerian methods as it is these methods that stand to benefit the most from our decomposi- tion approach. Nevertheless, complementary discussions surrounding the computation of the maximal reachable tube are provided in Section 4.5.5. In the new coordinates z = T−1x, T = T1T2, the subsystem dynamics are governed by ż1 = ( A11 −A12L− δMF (Z) ) z1 +B1u, (4.41) ż2 = ( A22 + LA12 + δF (Z)M ) z2 +B2u+ δF (Z)z1 (4.42) with δF (Z) = 0 when the standard Riccati transformation yields disjoint input, and B2 = 0 when the modified Riccati transformation is employed. Denote the two subspaces of Rn in which the subsystems evolve as S1 := Rk and S2 := Rn−k. (4.43) Algorithm 4.1 computes the (minimal) reachable tube in lower dimensions: 63 4.4. Reachability in Lower Dimensions Algorithm 4.1 Reachability in lower dimensions (Riccati-Based) 1: Zτ ← T−1K 2: Z iτ ← ProjSi(Zτ ), ∀i ∈ {1, 2} . project onto ith subspace For upper subsystem: 3: Z1[0,τ ] ← Reach[[0,τ ](Z1τ ,U) For lower subsystem: 4: Treat δF (Z)z1 as disturbance (existentially quantified) 5: ζ ← ‖z1‖ ≡ supv∈Z1 [0,τ ] ‖v‖ 6: Compute upper-bound ‖δF (Z)z1‖ ≤ ‖δF (Z)‖ζ 7: Z2[0,τ ] consrv.←− Reach][0,τ ](Z2τ ,B(‖δF (Z)‖ζ)) 8: return(Z1[0,τ ],Z2[0,τ ]) When the standard Riccati transformation is used to obtain the sub- systems, steps 4–7 of Algorithm 4.1 are simply replaced with Z2[0,τ ] ← Reach[[0,τ ](Z2τ ,U). Notice that if the error due to projections can be ignored, reachability in the upper subspace is “exact” in the sense that Z1[0,τ ] ≡ ProjS1(T−1Reach[[0,τ ](K,U)). This is generally not true in the lower subspace. The unidirectional cou- pling between the two subsystems is treated as a disturbance to the lower subsystem, resulting in a conservative reachability computation in that sub- space. The computed reachable tube in the lower subspace is a guaranteed over-approximation of the projection of the actual reachable tube in that subspace; Reach][0,τ ](Z2τ ,B(‖δF (Z)‖ζ)) ⊇ ProjS2(T−1Reach[[0,τ ](K,U)). The parameter δ = δ∗ is precomputed so as to minimize ‖δF (Z)‖. How- ever, when the reachability horizon in Step 3 of Algorithm 4.1 is large, the magnitude of the input to the lower subsystem (whose upper-bound is di- rectly proportional to supv∈Z1 [0,τ ] ‖v‖) may become so large as to warrant reachability in that subspace over sub-intervals of [0, τ ] in a similar fashion to [36]. For N := τ/q, N ∈ N time steps each of length q ∈ R+, we have Z2[0,τ ] = N−1⋃ i=0 Z2[iq, (i+1)q] (4.44) 64 4.4. Reachability in Lower Dimensions and, by the semi-group property, Z2[iq, (i+1)q] = Reach][0,q] (Z2[(i+1)q, (i+2)q],B(‖δF (Z)‖ζN−1−i)), (4.45) where Z2[τ,+∞) . = Z2τ and ζi is the supremum of the reachable tube in the upper subspace at intermediate time steps. This holds since the input to the lower subsystem is a disturbance and therefore all quantifiers in reachability analysis of this subsystem are existential. Recording ζi at each time step (rather than at the end of the reachability horizon) allows for a gradual incrementing of the disturbance to the lower subsystem. Thus, using the sequence {ζi}N−1i=0 and executing Algorithm 4.1 with sub-intervals according to (4.44)–(4.45), we can compute a less conservative reachable tube in the lower subspace. A guaranteed over-approximation of the actual reachable tube of the full- order system in X can be obtained using Lemma 3.1 as in the Schur-based case, i.e. via T ( (Z1[0,τ ] × S2) ∩ (S1 ×Z2[0,τ ]) ) ⊇ Reach[[0,τ ](K,U). (4.46) 4.4.1 Formulating an Upper-Bound on Conservatism of Z2[0,τ ] Consider the computed reachable tube Z2[0,τ ] := Reach][0,τ ](Z2τ ,B(‖δF (Z)‖ζ)) in the lower subspace. Take z̄2 ∈ Z2τ and suppose D[0,τ ] is the set of mea- surable functions from [0, τ ] to B(‖δF (Z)‖ζ). There exists an admissible input ϑ(·) ∈ D[0,τ ] such that in positive time using time-reversed dynamics we have z2 := e −τΩz̄2 − ∫ τ 0 e−(τ−r)Ωϑ(r)dr, (4.47) z2 ∈ Z2[0,τ ] (4.48) 65 4.4. Reachability in Lower Dimensions with Ω := A22 + LA12 + δF (Z)M . Therefore, as in [36], ‖z2 − e−τΩz̄2‖ ≤ ∫ τ 0 e(τ−r)‖Ω‖‖δF (Z)‖ζdr (4.49) = eτ‖Ω‖ − 1 ‖Ω‖ ‖δF (Z)‖ζ. (4.50) Notice that, eτ‖Ω‖ − 1 ‖Ω‖ ‖δF (Z)‖ζ = ( lim M→∞ M∑ i=1 τ i‖Ω‖i−1 i! ) ‖δF (Z)‖ζ (4.51) ≤ ( lim M→∞ M∑ i=1 τ i ( σ(Ω) √ nk )i−1 i! ) ‖δF (Z)‖ζ (4.52) =: η[0,τ ] (4.53) with σ(·) the largest singular value operator, and nk := n−k the dimension of the ETUC subsystem. Due to linearity and time-invariance of the dynamics, the (possibly non-convex) backward reachable tube can be expressed as Z2[0,τ ] ⊆ ⋃ t∈[0,τ ] e−ΩtZ2τ ⊕ B(η[0,τ ]). (4.54) The right-hand side of (4.54) provides an upper-bound on how much Z2[0,τ ] can grow in backward time in terms of the reachability horizon τ , the choice of nk, the magnitude of the unidirectional coupling ‖δF (Z)‖, the supremum of the reachable tube in the upper subspace ζ = supv∈Z1 [0,τ ] ‖v‖, and the largest singular value σ(Ω) of the lower subsystem. When reachability is performed over sub-intervals, a tighter upper-bound can be formulated by replacing the right-hand side of (4.54) with Z2[0,τ ] := N−1⋃ i=0 ⋃ t∈[iq,(i+1)q] e−ΩtZ2[(i+1)q,(i+2)q] ⊕ B(η[iq,(i+1)q]) . (4.55) As in [36], it can be shown that the quality of this upper-bound is good in 66 4.4. Reachability in Lower Dimensions the Hausdorff distance, in that limq→0 distH(Z2[0,τ ],Z 2 [0,τ ])→ 0. Determining whether Z2[0,τ ] itself is a close over-approximation to ProjS2(T−1Reach[[0,τ ](K,U)) is much more involved and remains an open problem. We expect that per- forming reachability over sub-intervals, choosing nk appropriately, and min- imizing the disturbance (uncertainty) magnitude as much as possible could reduce the conservatism considerably. 4.4.2 The Effect of Dimension on Magnitude of Uncertainty Since the worst-case unidirectional coupling ‖δF (Z)‖ contributes to the uncertainty (i.e. disturbance input) in the reachability computation for the lower subsystem, we examine the potential affect of the system dimension on a) the magnitude of the unidirectional coupling and b) the amount of time consumed by the decomposition process. Although entirely dependent on the particular system to which the mod- ified Riccati transformation is applied, the following empirical test can serve as a rough measure: For a given dimension n, we generated 10 randomized systems, applied the decomposition method to each n-dimensional system with nk = n 2 , and recorded ‖δF (Z)‖ and the computational time. We re- peated this for n = 4, 6, 8, 10, 16, as shown in Figure 4.2. While for higher dimensions, the computational time and the magnitude of the unidirectional coupling show an increasing trend in their average values, there is significant variance. In addition, the time required for the decomposition process, even for the highest dimension, is still negligible compared to the time required for actual reachability computations. 4.4.3 Conservatism Due to Projection We have assumed that the target set in the transformed coordinates is (or is close to) a direct product of the sets in the subsystems’ subspaces, meaning that the error due to projections of the target set onto the lower dimen- sional subspaces can be ignored. This assumption does not generally hold. If the target set in the new coordinate space is far from being axis-aligned, the projections contribute to the conservatism of the reachable tube over- 67 4.4. Reachability in Lower Dimensions 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 0 0.5 1 1.5 2 2.5 3 3.5 4 Decomposition time (s) M ag ni tu de o f u ni di re ct io na l c ou pl in g n=4 n=6 n=8 n=10 n=16 Figure 4.2: The worst-case unidirectional coupling and the time required to compute the transformed matrices for randomly generated systems of dimension n = 4, 6, 8, 10, 16 show significant variance. Average values are marked by ‘x’. approximation. A similar argument holds for the back-projection of the subsystem reachable tubes, since we attempt to reconstruct a higher di- mensional object from lower dimensional entities. Unfortunately, loss of information is inherent in any projection operation and, to a great extent, cannot be avoided. 4.4.4 Implications of Computing the Riccati-Based Reachable Tube Safety Verification Let I ⊂ X be the set of all possible initial states of (2.12). It follows from Definition 2.3 that, for a given unsafe target set K, the system is safe if and only if the backward minimal reachable tube does not intersect I: R[0,τ ] := Reach[[0,τ ](K,U) ∩ I = ∅. (4.56) Let R[0,τ ] and Z[0,τ ] denote the computed Riccati-based reachable tube in the original and transformed coordinates, respectively. Proposition 4.4. Z[0,τ ] ∩ T−1I = ∅ ⇐⇒ R[0,τ ] ∩ I = ∅ =⇒ R[0,τ ] ∩ I = ∅. 68 4.4. Reachability in Lower Dimensions Proof. ∀x ∈ R[0,τ ] ∀y ∈ I (R[0,τ ] ∩ I) = ∅ ⇐⇒ x 6= y ⇐⇒ T−1x 6= T−1y ⇐⇒ T−1R[0,τ ]∩T−1I = ∅. By Lemma 3.1 we have thatR[0,τ ] ⊆ R[0,τ ]. Thus, R[0,τ ] ∩ I = ∅ =⇒ R[0,τ ] ∩ I = ∅. Therefore, if the Riccati-based reachable tube does not intersect I in either coordinates, the system is safe. A simpler (but more conservative) sufficient condition for safety can be given as: Proposition 4.5. 2∧ i=1 ( Z i[0,τ ] ∩ ProjSi(T−1I) = ∅ ) =⇒ R[0,τ ] ∩ I = ∅. Proof. (Z1[0,τ ]∩ProjS1(T−1I) = ∅)∧(Z2[0,τ ]∩ProjS2(T−1I) = ∅) =⇒ (Z1[0,τ ]× S2)∩ (S1 ×Z2[0,τ ])∩ T−1I = ∅ ⇐⇒ Z[0,τ ] ∩ T−1I = ∅ =⇒ R[0,τ ] ∩ I = ∅. That is, if the computed Riccati-based reachable tubes for each subsys- tem in the new coordinates do not intersect the projections of the trans- formed I, then the system is safe in the original coordinates. This can be used to prove safety using the lower dimensional reachable tubes when reconstructing and storing the full-order state space (or reachable tube) is costly. Safety-Preserving Control Synthesis Given the unsafe target set K and the computed Riccati-based reachable tube R[0,τ ], let Rc[0,τ ] denote an under-approximation of the finite-time via- bility kernel V iab[0,τ ](Kc,U) under S. Then using the optimal control laws precomputed during the reachability analysis (e.g. via [90]), one can con- struct a feedback controller as in [81] on the boundaries of Rc[0,t], t ∈ [0, τ ], that keeps the trajectories of S within Kc (i.e. within safety) over the horizon [0, τ ]. Hence, although possibly conservative, computing the reachable tube through the Riccati-based approach can be a powerful tool to guarantee safety in safety-critical systems that have moderately high dimensions. 69 4.5. Numerical Examples 4.5 Numerical Examples We employ Level Set Toolbox (LS) [90] for the computation of the minimal reachable tubes (and the viability kernels). The use of our Riccati-based approach for maximal reachability analysis is discussed and an example is provided that employs Ellipsoidal Toolbox (ET) [73]. All computations are performed on a dual core Intel-based computer with 2.8 GHz CPU, 6 MB of L2 cache and 3 GB of RAM running single-threaded 32-bit Matlab 7.5. 4.5.1 Arbitrary 4D System Consider an LTI system ẋ = Ax+Bu with A = 1.5072 3.3984 0.1300 −0.0884 5.0644 −2.6683 0.0227 0.1689 0.1156 −0.1863 0.5686 0.2648 −0.0808 0.0229 0.4915 0.5949 , B = −0.7433 −2.2528 −0.9075 0.6036 and u ∈ R, |u| ≤ 0.1. We define a target (unsafe) set K such that in the transformed coordinate space Zτ = {z ∈ R4 | ‖z‖2 ≤ 0.2, z = T−1x, x ∈ K} where ‖·‖2 is the Euclidean norm and T is the transformation matrix obtained through the presented modified Riccati method. We decompose this system into two 2D subsystems. Sweeping through the real line, a nearly optimal δ∗ ≈ +50 that minimizes the worst-case uni- directional coupling between the subsystems is found in less than a second. Equations (4.27) and (4.37) converge to their fixed-point solutions in less than a dozen iterations. The system matrices in the new coordinate space are A′′ = 1.5362 3.4622 0 0 4.9377 −2.7030 0 0 −1.8075 −0.0193 0.5727 0.2612 1.2341 0.2918 0.4858 0.5964 , B′′ = −0.7433 −2.2528 0 0 . Reachability calculations are performed over a grid with 41 nodes in 70 4.5. Numerical Examples Figure 4.3: Series of 3D snapshots of the Riccati-based over-approximation (transparent light) vs. actual (solid dark) minimal reachable (unsafe) tube in the transformed coordinate space for Exam- ple 4.5.1. each dimension for τ = 3 s. The full-order reachable tube (10144.80 s com- putational time) is over-approximated by the Riccati-based reachable tube (3.98 s computational time, including calculation of δ∗, transformation ma- trices, the decomposition, and projections) as shown in Figure 4.3. The loss of accuracy (due to treating the unidirectional coupling as dis- turbance to the lower subsystems and assembling the full-order reachable tube from projections) can be quantified in terms of the Hausdorff distance between the full-order and the Riccati-based reachable tubes: Since a set in LS is described by a signed distance function whose magnitude at any point in the state space is the minimum distance to the boundary of that set, it is straightforward to compute the Hausdorff distance between any given two sets. For this example the Hausdorff distance is 0.18. In addition, a volumetric measure of the inaccuracy can also be quantified in terms of the 71 4.5. Numerical Examples ratio between the number of grid points contained in the difference of the two sets and the number of grid points in the Riccati-based reachable tube. This ratio is found to be 38.4%. 4.5.2 4D Cart with Two Inverted Pendulums Consider the linearized model of a cart with two separately mounted inverted pendulums from [51, Ex. 2.2.1] with l1 = 30, l2 = 35. The state vector x ∈ R4 consists of angular displacement of each inverted pendulum from vertical and the corresponding angular velocities; the control input u ∈ R arises from a force applied to the cart such that |u| ≤ 10. The system matrices are A = 0 1 0 0 0.3920 0 −0.0327 0 0 0 0 1 0.0560 0 0.2753 0 , B = 0 −0.0033 0 −0.0005 . We choose a non-convex target (unsafe) set K such that in the trans- formed coordinates we have Zτ = {z ∈ R4 | ‖z‖ ≥ 0.5, z = T−1x, x ∈ K}. We seek to identify the set of states for which there exists a bounded con- trol law that keeps the system trajectories contained in Zcτ . The safety- preserving control synthesized through LS provides a guarantee that the pen- dulums’ angular displacement will not exceed an infinity norm ball around their upright positions, despite control saturation. Similar “envelope protec- tion” problems arise in other domains, including aircraft flight management systems and anesthesia automation, among others. We decompose this system into two 2D subsystems, with the unidirec- tional coupling determined by the solution L = −B2B†1 regardless of the 72 4.5. Numerical Examples Figure 4.4: Riccati-based (solid dark) vs. actual (transparent light) viability kernels in the transformed coordinate space for Example 4.5.2. The minimal reachable tube and its over-approximation are the non-convex complements of these objects. value of δ. The system matrices in the new coordinate space become A′′ = 0 0.9524 0 0 0.3920 0 0 0 0 0.1429 0 1.0500 0 0 0.2800 0 , B′′ = 0 −0.0033 0 0 . Reachability calculations are performed over a grid with 41 nodes in each dimension for τ = 3 s. Figure 4.4 shows the computed viability kernels (safe) as the area inside the shaded regions. The computation time for the actual and the transformation-based reachable tubes were 1098.48 s and 4.27 s, respectively. The Hausdorff distance between these two sets is 0.21 and the Riccati-based viability kernel covers 74% of the volume of the full- order set. 73 4.5. Numerical Examples 4.5.3 Arbitrary 6D System Consider the system ẋ = Ax+Bu with A = 3.3155 0.7768 2.4455 0.0028 −0.0094 −0.0097 0.6320 −1.4796 −2.3001 −0.0370 0.0322 −0.0112 −0.1047 −0.3522 0.6578 0.0282 −0.0621 −0.0214 0.0344 0.0360 0.0091 0.1885 0.0518 0.3567 −0.0140 −0.0225 −0.0012 0.2746 0.0866 0.1567 −0.0106 −0.0215 −0.0082 0.0814 0.0887 0.1182 , B = [ 0.1469 −0.7988 −2.3854 −0.0054 1.3260 −0.1623 0.2657 2.4582 −0.3955 −0.1403 −0.1187 0.3601 ]T and u ∈ R2, ‖u‖ ≤ 0.05. We decompose this system into two 3D sub- systems using the modified Riccati transformation with δ∗ ≈ −199 (see Appendix B.3). A non-convex target set K is chosen such that in the trans- formed coordinates Z iτ = ⋃3 j=1 Pj , ∀i ∈ {1, 2}, where Pj = {zi ∈ R3 | Czi − bj ≤ 0} with CT = 1 −1 0 0 0 00 0 1 −1 0 0 0 0 0 0 1 −1 , bT1 = [ 0.5 0.3 0.1 0.1 −0.2 0.4 ] , bT2 = [ −0.1 0.3 0.1 0.1 0.4 0.4 ] , bT3 = [ 0.5 0.3 0.1 0.1 0.4 −0.2 ] , and z = [zT1 , z T 2 ] T = T−1x, x ∈ K. Reachability computations are performed over a grid with 71 nodes in each dimension for τ = 2 s, as shown in Figure 4.5. The full-order 6D (minimal) reachable tube is the intersection of the back-projection of these 3D reachable tubes. The overall computation time was 489.89 s, whereas the actual reachable tube is prohibitively computationally expensive to compute with LS for any meaningful grid resolution. In addition, only 28 MB of RAM was used in the Riccati-based calculations, whereas computation of the full- order reachable tube would require over 4 TB of RAM (well beyond the capabilities of today’s technology). 74 4.5. Numerical Examples (a) Upper subspace (b) Lower subspace Figure 4.5: The Riccati-based reachable tube (unsafe) in the transformed coordinates for Example 4.5.3. 4.5.4 Comparison With Schur-Based Decomposition (Chapter 3) In Chapter 3 we presented a Schur-based decomposition technique that is applicable to almost any LTI system (subject to a mild assumption on the input-to-state map). In contrast, the decomposition method presented here is based on two nonsymmetric algebraic equations. The existence of so- lutions to these algebraic equations, however, is limited by a number of conditions on system matrices and is therefore heavily problem dependent. Indeed, as pointed out earlier, the conditions are more likely to be satisfied as the ill-conditioning of the original system matrices increases—e.g., for two-time-scale systems. (Figure 4.6 shows the fraction of tests on randomly generated systems for which a solution existed.) However, when the alge- braic Riccati equations converge, the resulting subsystems could potentially yield less conservative reachable tube over-approximations than in the case of the Schur-based decomposition. Consider a simple constrained 2D system with A = [−2.0228 0.9732 −0.3695 0.0893 ] , B = [ 0.96000.1372 ]. Applying the modified Riccati transformation results in A′′ric = [−1.8360 0 −0.0875 −0.0975 ] , B′′ric = [ 0.96000 ]. Notice that the “fast” eigenvalue λ1 = −1.8360 is assigned to the controlled subsystem and that the mag- 75 4.5. Numerical Examples 4 6 8 10 16 0 50 100 Dimension Su cc es s ra te [% ] 93.38% 88.08% 84.11% 79.47% 68.87% Figure 4.6: Fraction of a randomized test for which a Riccati solution existed. Success rate is shown as a percentage of 150 ran- domly generated two-time-scale systems for each dimension n = 4, 6, 8, 10, 16 with nk = n 2 . The A matrix entries of each system are drawn from a normal distribution with mean 0 and standard deviations 2 and 0.2 for A11 and A22 blocks respec- tively. nitude of the unidirectional coupling is ‖δF (Z)‖ = 0.0875. To retain the same eigen-structure (that is, the controlled subsystem be associated with the fast eigenvalue and the uncontrolled and perturbed subsystem with the slow eigenvalue λ2 = −0.0975) using the Schur-based decomposition yields A′′sch = [−0.0975 −0.1276 0 −1.8360 ] , B′′sch = [ 0−0.7948 ] . The corresponding subsystems are subject to a unidirectional coupling that is 46% larger in magnitude than in the Riccati-based case. Suppose the target set is the Euclidean unit disk. This set is reshaped similarly under both transformations (Figure 4.7). Since the magnitude of the input-to-state map in the Schur-based case is smaller than that in the Riccati-based case, the optimal Hamiltonian for the reachability analysis of its controlled subsystem is also smaller. Therefore the reachable tube of the controlled subsystem in the Schur-based case is larger than in the Riccati-based case. In both decomposition techniques, the supremum of the reachable tube of the controlled subsystem as well as the unidirectional coupling between the subsystems are treated as disturbance in reachability computation of the ETUC subsystems. Since this disturbance has a larger magnitude in the Schur-based case, the overall reachable tube computation is more conservative than in the Riccati-based case. Therefore the Riccati- based decomposition is the superior approach in this example. 76 4.5. Numerical Examples −1.5 −1 −0.5 0 0.5 1 1.5−1.5 −1 −0.5 0 0.5 1 1.5 Schur Riccati Figure 4.7: The unit disc target set under Riccati-based and Schur-based transformations. The corresponding reachable tube using the Riccati-based approach will be less conservative. One should note, however, that no universal conclusions can be drawn from this particular example: In general, it is the problem under study (i.e. the system matrices, the shape of the target set in the transformed coordinate space, the effect of projections, etc.) that determines which de- composition method is the better choice. 4.5.5 The Decomposition and Maximal Reachability Analysis With an unsafe target set, the Riccati-based approach presented here can also be used (in conjunction with both Eulerian and Lagrangian techniques) to over-approximate the maximal reachable tube, where the input u to the original system is considered as “disturbance” or “uncertainty” and is exis- tentially quantified. This is done by replacing Step 3 in Algorithm 4.1 with Z1[0,τ ] ← Reach][0,τ ](Z1τ ,U) (and Steps 4–7 with Z2[0,τ ] ← Reach][0,τ ](Z2τ ,U) if the standard Riccati transformation is used to obtain the subsystems). As expected, we have ProjS1(T −1Reach][0,τ ](K,U)) ≡ Reach][0,τ ](Z1τ ,U), (4.57) and ProjS2(T −1Reach][0,τ ](K,U)) ⊆ Reach][0,τ ](Z2τ ,B(‖δF (Z)‖ζ)). (4.58) 77 4.5. Numerical Examples The over-approximation of the actual maximal reachable tube of the full- order system can be obtained similarly to Lemma 3.1 as T ( (Z1[0,τ ] × S2) ∩ (S1 ×Z2[0,τ ]) ) ⊇ Reach][0,τ ](K,U). (4.59) While Algorithm 4.1 with the above modifications ensures that the full- order (unsafe) maximal reachable tube is over-approximated, the results may be too conservative. Recall that the modified Riccati transformation is best suited to systems that are two-time-scale and ill-conditioned. The fast eigenvalues are assigned to the upper subsystem whereas the slow eigenval- ues to the lower subsystem. Such an eigenvalue allocation is advantageous when computing the minimal reachable tube. In maximal reachability anal- ysis, however, since the input is existentially quantified, the reachable tube in the upper subspace can grow significantly. As a result, the input to the lower subsystem, whose upper bound is directly proportional to the supre- mum of the reachable tube in the upper subspace, may become excessively large. Consequently, the reachable tube in the lower subspace may be overly conservative. Note that a conservative over-approximation may be justified when the constraints are severely non-convex and/or cannot be adequately approx- imated by compact convex shapes (or a union of a few such sets), and therefore the computationally intensive Eulerian methods must be used to compute the maximal reachable tube. With convex constraints, however, the Lagrangian methods can be utilized to compute the full-order maximal reachable tube with great accuracy and efficiency. This will also circumvent the errors introduced via projections, etc. as compared to the case in which the Riccati-based approach is employed. Consider a 4D example ẋ = Ax+Bu, u ∈ R, |u| ≤ 0.1, with A = 0.5990 −0.2333 −0.0244 0.0121 −0.4553 −0.2107 −0.0076 −0.0041 −0.0091 0.0128 −0.1749 −0.1768 0.0309 −0.0229 −0.0105 −0.0777 , B = 0.6846 2.4813 0.3583 −1.5195 . 78 4.6. Summary and Further Discussions −2 −1 0 1 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 z1 z 2 (a) (z1, z2) subspace −2 −1 0 1 2 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 z3 z 4 (b) (z3, z4) subspace Figure 4.8: Projections of the full-order maximal reachable tube (dark/blue) vs. the Riccati-based maximal reachable tube (light/red) in the transformed coordinates, computed using ET. (The two sets are indistinguishable in (z1, z2) subspace.) The target set in the transformed coordinates is the Euclidean unit ball. We use the Riccati-based method to decompose this system into two 2D subsystems (see Appendix B.3). The computation time for this decompo- sition (including the search for δ∗ ≈ 12.6) was 0.42 s. To over-approximate the maximal reachable tube of this system we employ ET [73] and perform lower dimensional reachability over sub-intervals for τ = 1 s. Figure 4.8 com- pares the projections of the full-order reachable tube (computed in 1.08 s) and the lower-dimensional reachable tubes (collectively computed in 4.1 s) in the transformed coordinates. The greater computational time for the lower- dimensional reachability may be associated with the for-loops and other demanding computations involved in a sub-interval calculation. 4.6 Summary and Further Discussions In this chapter we presented our second decomposition method for reach- ability analysis of LTI systems based on the Riccati transformation. This decomposition has considerable potential for reducing the computational complexity in reachability calculations, particularly for Eulerian methods. 79 4.6. Summary and Further Discussions For the case in which the input was non-disjoint across resulting subsys- tems, a modified Riccati transformation was proposed. The extension of this transformation-based reachability approach to switched/hybrid systems with LTI continuous dynamics can be easily achieved following the procedure described in Section 3.5. Severely non-convex constraints and/or the need to compute the mini- mal reachable tube and the viability kernel as well as their corresponding safety-preserving control laws warrant the use of computationally intensive Eulerian methods. Despite inevitable conservatism e.g. due to the use of projection, our approach provides a means to compute the minimal reach- able tube for higher dimensional systems for which these computationally intensive reachability tools were previously not applicable. It is possible (although uncommon) that the transformation matrix can become poorly-conditioned due to pseudoinverses and numerical algorithms involved, resulting in the target set in the transformed coordinates becom- ing too severely distorted under the linear map to be of any practical use. An upper-bound on the condition number in terms of the system matrices and the free parameter δ is provided in Appendix B.2. We are currently investigating possible remedies that would ensure a well-conditioned trans- formation matrix. Finally, it is worth emphasizing that the proposed modified Riccati trans- formation has application beyond reachability analysis. The technique can be used in any scenario in which decoupling of the dynamics as well as the input-to-state map is required. 80 Chapter 5 Set-Theoretic Methods: Lagrangian Algorithms for Viability1 In the following two chapters we will address Problem 2.2 and will present our second approach, based on set-theoretic methods, to reduce the com- plexity of the computation of the minimal reachable tube or the viability kernel. We will make use of the definitions and terminologies introduced in Section 2.5. Consider a given constraint/target set K ⊂ X . While our main focus in Chapters 3 and 4 was on approximating the minimal reachable tube Reach[[0,τ ](K,U) when K was deemed unsafe, in this chapter we will turn most of our attention to the case in which K is deemed safe and we seek to approximate the viability kernel V iab[0,τ ](K,U) of this set. As mentioned before, approximating the minimal reachable tube and the viability kernel are not mutually exclusive as these constructs are duals of one another. For example, a method that facilitates an under-approximation of the viability kernel of K automatically provides a means for the over-approximation of the minimal reachable tube for Kc. Our approach is based on drawing a connection between the viability kernel and maximal reachable sets. As mentioned extensively in previous chapters, the Eulerian schemes normally used to compute the viability kernel suffer from a complexity that is exponential in the dimension of the states. In contrast, the efficient and scalable Lagrangian methods compute maximal 1A version of this chapter has been published in [55]. 81 5.1. The Viability Kernel in Terms of Maximal Reach Sets reachable sets. We will show that under certain conditions these methods can be employed to conservatively approximate the viability kernel for pos- sibly high-dimensional systems. Significant reduction in the computational costs can be achieved since instead of a single calculation with exponential complexity one can perform a series of calculations with polynomial com- plexity. 5.1 Connection Between the Viability Kernel and Maximal Reachable Sets We begin the analysis by considering the continuous-time case first and then proceed to the discrete-time case. 5.1.1 Continuous-Time Systems Consider the case in which (2.1) is the continuous-time nonlinear system ẋ(t) = f(x(t), u(t)), x(0) = x0, t ∈ R+. (5.1) We will show that we can approximate V iab[0,τ ](K,U) using a nested se- quence of sets that are reachable in small sub-intervals of [0, τ ]. Computing an Under-Approximation of the Viability Kernel Assume that the vector field f is bounded by M in the norm ‖·‖. Consider a partition P ∈P([0, τ ]). We begin by defining an under-approximation of the state constraint set (Figure 5.1(a)): K↓(P ) := {x ∈ K | dist(x,Kc) ≥M‖P‖}. (5.2) We under-approximate K by a distance M‖P‖ so as to only consider the system’s state at discrete times t0, t1, . . . , tn. At a time t in the interval 82 5.1. The Viability Kernel in Terms of Maximal Reach Sets [ti, ti+1], a trajectory x(·) can travel a distance of at most ‖x(t)− x(ti)‖ ≤ ∫ t ti ‖ẋ(τ)‖dτ ≤M(t− ti) ≤M‖P‖ (5.3) from its initial location x(ti). As we shall see, formulating the subset (5.2) will ensure that the state does not leave K at any time during [0, τ ]. The set K↓(P ) defines the first step of our recursion. We then define a sequence of |P | sets recursively: K|P |(P ) = K↓(P ), (5.4a) Kk−1(P ) = K↓(P ) ∩Reach]tk−tk−1(Kk(P ),U) for k ∈ {1, . . . , |P |}. (5.4b) At each time step, we calculate the set of states from which you can reach Kk(P ), then intersect this set with the set of safe states (see Figure 5.1). The final set K0(P ) is an approximation of V iab[0,τ ](K,U). Note that the resulting set depends on our choice of a partition P of the time interval [0, τ ]. We claim that for any partition P , K0(P ) is an under-approximation. Proposition 5.1. Suppose that the vector field f : X × U → X is bounded on a set K ⊆ X . Then for any partition P of [0, τ ] the final set K0(P ) defined by the recurrence relation (5.4) satisfies K0(P ) ⊆ V iab[0,τ ](K,U). (5.5) 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, fix a partition P of [0, τ ] and take a point x0 ∈ K0(P ). By the construction of K0(P ), this means that for each k = 1, . . . , |P | there is some point xk ∈ Kk(P ) and an input uk : [0, tk−tk−1]→ U such that xk can be reached from xk−1 at time tk−tk−1 using input uk. Thus, taking the concatenation of the inputs uk, we get an input u : [0, τ ]→ U such that the solution x : [0, τ ]→ X to the initial value problem ẋ = f(x, u), x(0) = x0, satisfies x(tk) = xk ∈ Kk(P ) ⊆ {x ∈ K | 83 5.1. The Viability Kernel in Terms of Maximal Reach Sets (a) We define the initial under-approximation of the safe set K|P |(P )=K↓(P ) (b) We calculate the set of backward reachable states from K|P |(P ) (c) We intersect the backward reachable set with the initial set to get K|P |−1(P ) (d) Next, we calculate the set of backward reachable states from K|P |−1(P ) (e) Again, we intersect the backward reachable set with the initial set to get a new set K|P |−2(P ) (f) By repeating this pro- cess, we reach an under- approximation K0(P ) of the viability kernel. Figure 5.1: Iteratively constructing an under-approximation of V iab[0,τ ](K,U). 84 5.1. The Viability Kernel in Terms of Maximal Reach Sets dist(x,Kc) ≥ M‖P‖}. We claim that this guarantees that x(t) ∈ K for all t ∈ [0, τ ]. Indeed, any t ∈ [0, τ) lies in some interval [tk, tk+1). Since f is bounded by M , we have ‖x(t)− x(tk)‖ ≤M(t− tk) < M(tk+1 − tk) ≤M‖P‖. (5.6) Further, x(tk) ∈ Kk(P ) implies that dist(x(tk),Kc) ≥ M‖P‖. Combining these, we see that dist(x(t),Kc) ≥ dist(x(tk),Kc)− ‖x(t)− x(tk)‖ > M‖P‖ −M‖P‖ = 0 (5.7) and hence x(t) ∈ K. Thus, x0 ∈ V iab[0,τ ](K,U). Remark 5.1. For a large enough horizon τ , if Kk−1(P ) ≡ Kk(P ) for some k ∈ {1, . . . , |P |} with tk ∈ (0, τ ], then K0(P ) = Kk(P ) approximates the infinite-horizon viability kernel V iabR+(K,U). This set is also known as the maximal controlled-invariant subset [12] (see Section 2.1). Precision of the Approximation The approximation can be made to be arbitrarily precise by choosing a sufficiently fine partition. This is true in the sense that the union of the approximating sets K0(P ) taken over all possible partitions P of [0, τ ] is bounded between the viability kernels of K and its interior ◦K. Proposition 5.2. Suppose that the vector field f : X × U → X is bounded on a set K ⊆ X . Then we have V iab[0,τ ]( ◦K,U) ⊆ ⋃ P∈P([0,τ ]) K0(P ) ⊆ V iab[0,τ ](K,U). (5.8) In particular, when K is open,⋃ P∈P([0,τ ]) K0(P ) = V iab[0,τ ](K,U). (5.9) 85 5.1. The Viability Kernel in Terms of Maximal Reach Sets Kc ∂K x(·) dM‖P‖ K̊ x0 x(τ) Figure 5.2: Distance d > 0 from the boundary set ∂K. Partition P can be chosen such that M‖P‖ < d. Proof. The second inclusion in (5.8) follows directly from Proposition 5.1. To prove the first inclusion, take a state x0 ∈ V iab[0,τ ]( ◦K,U). 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, τ ]→ X is continuous so the function t 7→ dist(x(t),Kc) is continuous on the compact set [0, τ ]. Thus, we can define d > 0 to be its minimum value. Now take a partition P of [0, τ ] such that M‖P‖ < d (see Figure 5.2). We need to show that x0 ∈ K0(P ). First note that our partition P is chosen such that dist(x(t),Kc) > M‖P‖ for all t ∈ [0, τ ]. Hence x(tk) ∈ K|P |(P ) for all k = 0, . . . , |P |. To show that x(tk−1) ∈ Reach]tk−tk−1(Kk(P ),U) for all k = 1, . . . , |P |, consider the tokenization2 {uk}k of the input u corresponding to P . It is easy to verify that for all k, we can reach x(tk) from x(tk−1) at time tk − tk−1 using input uk. Thus, in particular, we have x0 = x(t0) ∈ Reach]t1−t0(K1(P ),U). So x0 ∈ K0(P ). Hence V iab[0,τ ]( ◦K,U) ⊆ ⋃P∈P([0,τ ])K0(P ). 2See Definition 2.11. 86 5.1. The Viability Kernel in Terms of Maximal Reach Sets 5.1.2 Discrete-Time Systems Consider the case in which (2.1) is the discrete-time nonlinear system x(t+ 1) = f(x(t), u(t)), x(0) = x0, t ∈ Z+. (5.10) Computing V iab[0,τ ](K,U) under this system is a particular case of the re- sults presented in Section 5.1.1. Define a sequence of sets recursively as Kn = K, (5.11a) Kk−1 = K ∩Reach]1(Kk,U) for k ∈ {1, . . . , n}, (5.11b) where τ = n and Reach]1(·, ·) is the unit time-step maximal reachable set. Proposition 5.3. Let K0 be the final set obtained from the recurrence re- lation (5.11). Then, V iab[0,τ ](K,U) = K0. (5.12) Proof. Without loss of generality we assume that the time variable t is in- teger valued. As a result, the tokenization of the input signal u is a discrete sequence {uk}k with uk := u(t) with t = k − 1 for k = 1, . . . , n. To show K0 ⊆ V iab[0,τ ](K,U), via recursion (5.11) we have that at each step k there exists uk such that xk−1 ∈ Kk−1 reaches xk ∈ Kk. Thus, x0 ∈ K0 implies there exists a concatenation u(·) = {uk}k ∈ U[0,τ ] such that x(t) ∈ K for all t ∈ [0, τ ]. Therefore, x0 ∈ V iab[0,τ ](K,U). To show V iab[0,τ ](K,U) ⊆ K0, take x0 ∈ V iab[0,τ ](K,U). There exists u(·) = {uk}k such that x(t) ∈ K for every t. Using the tokenization of {uk}k we can verify that for some uk we can reach xk := x(t+1) from xk−1 := x(t). Hence, xk−1 ∈ Reach]1(Kk,U) for all k ∈ {1, . . . , n}. In particular, for k = 1 we have x0 := x(0) ∈ Reach]1(K1,U). Thus, x0 ∈ K ∩ Reach]1(K1,U) = K0. Remark 5.2. Note that the above iterative scheme is closely related to the set-valued description of the discrete viability kernel presented in [18, 106] 87 5.2. Computational Algorithms and the recursive construction of the controlled-invariant subset for discrete- time systems presented in [11, 12, 61, 120]. 5.2 Computational Algorithms Any technique that is capable of computing the maximal reachable set can now be used to compute the viability kernel. Most currently avail- able Lagrangian methods yield an (under- and/or over-) approximation of the maximal reachable set. The viability kernel should not be over- approximated since an over-approximation would contain initial states for which the viability of the system is inevitably at stake. Thus, to correctly compute V iab[0,τ ](K,U) all approximations must be in the form of under- approximations. Every step of the recursions (5.4) and (5.11) involves a reachability com- putation and an intersection operation. Ideally, the sets that are being intersected should be drawn from classes of shapes that are closed under such an operation, e.g. polytopes. However, the currently available reacha- bility techniques that are based on polytopes (e.g. [76]) do not, in general, scale well with the dimension of the state. Moreover, the scalable reacha- bility techniques, such as the methods of zonotopes [40], ellipsoids [70, 74], and support functions [39], generate sets that may prove to be difficult to transform into a (under-approximating) polytope. For instance, one may compute a polytopic under-approximation of the reachable sets using their support functions based on the approach presented in [77]. However, that approach requires calculation of the facet representation of the resulting polytopes from their vertices before each intersection operation, which is known to be computationally demanding in higher dimensions. Recently, the authors in [54] introduced an efficient polytopic technique with a fixed complexity called bounded vertex representation that is capable of computing under-approximations of the maximal reachable set and the intersection of polytopes. Unfortunately, this method assumes that the dy- namics are affine (of the form Ax+b) which circumvents the difficulties that arise from linear transformations and Minkowski sums in the reachability 88 5.2. Computational Algorithms operation. 5.2.1 A Piecewise Ellipsoidal Approach Here we showcase our results using an efficient algorithm, based on the ellipsoidal techniques [70] implemented in Ellipsoidal Toolbox (ET) [73], that sacrifices accuracy in exchange for scalability. We consider the LTI system L(x(t)) = Ax(t) +Bu(t) (5.13) with A ∈ Rn×n and B ∈ Rn×m. As in (2.1), depending on whether the system evolves in continuous time (t ∈ R+) or discrete time (t ∈ Z+), L(·) denotes the derivative operator or the unit forward shift operator, respec- tively. An ellipsoid in Rn is defined as E(q,Q) := {x ∈ Rn | 〈(x− q), Q−1(x− q)〉 ≤ 1} (5.14) with center q ∈ Rn and shape matrix Rn×n 3 Q = QT 0. A piecewise ellipsoidal set is the union of a finite number of ellipsoids. Among many advantages, ellipsoidal techniques [70, 73] allow for an ef- ficient computation of under-approximations of the maximal reachable sets, making them a particularly attractive choice for the reachability computa- tions involved in our formulation of the viability kernel. Suppose K and U are (or can be closely under-approximated as) compact ellipsoids with nonempty interior. Consider the continuous-time case and the recursion (5.4). (The arguments in the discrete-time case are similar.) Given a partition P and some k ∈ {1, . . . , |P |}, let Kk(P ) = E(xδ, Xδ) ⊂ X . As in [71], with N := {v ∈ Rn | 〈v, v〉 = 1} and δ := tk − tk−1 we have Reach]δ−t(Kk(P ),U) = ⋃ `δ∈N E(xc(t), X−` (t)), ∀t ∈ [0, δ], (5.15) where xc(t) and X−` (t) are the center and the shape matrix of the internal approximating ellipsoid at time t that is tangent to Reach]δ−t(Kk(P ),U) in 89 5.2. Computational Algorithms the direction `(t) ∈ Rn. For a fixed `(δ) = `δ ∈ N , the direction `(t) is obtained from the adjoint equation ˙̀(t) = −AT`(t). The center xc(t) (with xc(δ) = xδ) and the shape matrix X − ` (t) (with X − ` (δ) = Xδ) are determined from differential equations described in [72]. (cf. [74] for their discrete-time counterparts.) In practice, only a finite number of directions is used for the maximal reachable set computations. Let M be a finite subset of N . Then, Reach]δ−t(Kk(P ),U) ⊇ ⋃ `δ∈M E(xc(t), X−` (t)), ∀t ∈ [0, δ]. (5.16) Note that the under-approximation in (5.16) is in general an arbitrarily shaped, non-convex set. Performing our desired operations on this set while maintaining efficiency may be difficult, if not impossible. Instead consider the final backward reachable set Reach]δ(Kk(P ),U) and let Reach ][`δ] δ (Kk(P ),U) denote the maximal reachable set corresponding to a single terminal direction `δ := `(δ) ∈M. We have that Reach ][`δ] δ (Kk(P ),U) = E(xc(0), X−` (0)) ⊆ ⋃ `δ∈M E(xc(0), X−` (0)) ⊆ Reach]δ(Kk(P ),U). (5.17) Therefore, the reachable set computed for a single direction is an ellipsoidal subset of the actual reachable set. Let # : 2X → 2X denote a set-valued function that maps a set to its maximum volume inscribed ellipsoid. Algorithms 5.1 and 5.2 compute a piecewise ellipsoidal under-approximation of V iab[0,τ ](K,U) for continuous- time and discrete-time systems, respectively. 90 5.2. Computational Algorithms Algorithm 5.1 Piecewise ellipsoidal approximation of V iab[0,τ ](K,U) (continuous-time) 1: Choose P ∈P([0, τ ]) . Affects precision of approximation. 2: K|P |(P )← K B(M‖P‖) . Find {x ∈ K | dist(x,Kc) ≥M‖P‖}. 3: K∗0 (P )← ∅ 4: while M 6= ∅ do 5: l← `τ ∈M 6: k ← |P | 7: while k 6= 0 do 8: if Kk(P ) = ∅ then 9: K0(P )← ∅ 10: break 11: end if 12: G ← Reach][l]tk−tk−1(Kk(P ),U) . Compute the ellipsoidal under-approximation of the max- imal reach set along the direction l. 13: Kk−1(P )← #(K|P |(P ) ∩ G) . Find the max volume inscribed ellipsoid in K|P |(P ) ∩ G. 14: k ← k − 1 15: end while 16: K∗0 (P )← K∗0 (P ) ∪K0(P ) 17: M←M\{l} 18: end while 19: return (K∗0 (P )) 91 5.2. Computational Algorithms Algorithm 5.2 Piecewise ellipsoidal approximation of V iab[0,τ ](K,U) (discrete-time) 1: Kn ← K 2: K∗0 ← ∅ 3: while M 6= ∅ do 4: l← `τ ∈M 5: k ← n 6: while k 6= 0 do 7: if Kk = ∅ then 8: K0 ← ∅ 9: break 10: end if 11: G ← Reach][l]1 (Kk,U) 12: Kk−1 ← #(Kn ∩ G) 13: k ← k − 1 14: end while 15: K∗0 ← K∗0 ∪K0 16: M←M\{l} 17: end while 18: return (K∗0 ) Proposition 5.4. For a given partition P ∈P([0, τ ]), let K∗0 (P ) be the set generated by Algorithm 5.1. Then, K∗0 (P ) ⊆ V iab[0,τ ](K,U). (5.18) Proof. Let K̃0(P ) denote the final set constructed recursively by (5.4). Also, for a fixed direction l, let K [l] 0 (P ) denote the set produced at the end of each outer loop in Algorithm 5.1. Notice that via (5.17), for every l ∈ M, K [l] 0 (P ) ⊆ K̃0(P ). Therefore, ⋃ l∈MK [l] 0 (P ) ⊆ K̃0(P ). Thus, K∗0 (P ) =⋃ l∈MK [l] 0 (P ) ⊆ V iab[0,τ ](K,U). Remark 5.3. A similar argument holds for the discrete-time case in Algo- rithm 5.2, i.e. K∗0 ⊆ V iab[0,τ ](K,U). 92 5.2. Computational Algorithms #(·): Computing the Maximum Volume Inscribed Ellipsoid Notice that in the continuous-time case, the sets Y := K|P |(P ) and G := Reach ][l] tk−tk−1(Kk(P ),U) are compact ellipsoids for every l ∈M, P ∈P([0, τ ]), and k ∈ {1, . . . , |P |}. Similarly in the discrete-time case, Y := Kn and G := Reach ][l] 1 (Kk,U) are compact ellipsoids for every l ∈M and k ∈ {1, . . . , n}. Their intersection is, in general, not an ellipsoid but can be easily under- approximated by one. The operation #(·) under-approximates this intersec- tion by computing the maximum volume inscribed ellipsoid in Y ∩ G. The result is an ellipsoid that, while aiming to minimize the accuracy loss, can be used directly as the target set for the reachability computation in the subsequent time step. Let us rewrite the general ellipsoid as E(q,Q) = {Hx + q | ‖x‖2 ≤ 1} with H = Q 12 . Assume Y ∩ G 6= ∅ and suppose Y = E(q1, Q1) and G = E(q2, Q2). Following [15], the computation of the maximum volume inscribed ellipsoid in Y ∩ G (a readily-available feature in ET) can be cast as a convex semidefinite program (SDP): minimize H∈Rn×n, q∈Rn, λi∈R log detH−1 (5.19a) subject to λi > 0 (5.19b)1− λi 0 (q − qi) T 0 λiI H q − qi H Qi 0, i = 1, 2. (5.19c) Using the optimal values for H and q, we will have #(Y ∩ G) = E(q,HTH). (5.20) Loss of Accuracy A set generated by Algorithms 5.1 or 5.2 could be an inaccurate approxima- tion of V iab[0,τ ](K,U), especially for large time horizons. The loss of accu- racy is mainly attributed to the function #(·), the under-approximation of the intersection at every iteration with its maximum volume inscribed ellip- 93 5.2. Computational Algorithms soid. This approximation error propagates through the algorithms making them subject to the “wrapping effect” (cf. [66]). In the continuous-time case, the quality of approximation is also affected by the choice of time interval partition (Proposition 5.2). Choosing a finer partition increases the quality of approximation. However, doing so would also require a larger number of intersections to be performed in the inter- mediate steps of the recursion. As such, one would expect that the error generated by #(·) would be amplified. Luckily, since with a finer partition the reachable sets change very little from one time step to the next, the intersection error at every iteration becomes smaller. The end result is a smaller accumulated error and therefore a better approximation, at least experimentally. We show this using a trivial example: Consider the double integrator ẋ(t) = [ 0 1 0 0 ] x(t) + [ 0 1 ] u(t) (5.21) subject to ellipsoidal constraints u(t) ∈ U := [−0.25, 0.25] and x(t) ∈ K := E(0, [ 0.25 00 0.25 ]), ∀t ∈ [0, 1]. We employ eight different partitions P of the time interval such that we have |P | = 13, 21, 34, 55, 89, 144, 233, 377 and ‖P‖ = 1/|P |. The linear vector field is bounded on K in the infinity norm by M = ‖[ 0 10 0 ]‖ supx∈K‖x‖+ ‖[ 01 ]‖ supu∈U‖u‖ = 0.75. Thus, in Algorithm 5.1, K|P |(P ) = K B(0.75×‖P‖). A piecewise ellipsoidal under-approximation of V iab[0,1](K,U) for every partition P (with |M| = 10 randomly chosen initial directions) is shown in Figure 5.3. Notice that as |P | increases, the fidelity of approximation improves. A plot of the error in the under-approximation as a function of |P | is provided in Figure 5.4. Forming the Under-Approximation K↓(P ) in the Continuous-Time Case While a straightforward method to construct the under-approximationK↓(P ) of the set K in Algorithm 5.1 for a fixed P ∈ P([0, τ ]) is to erode K by a ball of radius M‖P‖ (for a given uniform bound M on f(x, u) = Ax+Bu), 94 5.2. Computational Algorithms -0.5 0 0.5 -0.5 0 0.5 -0.5 0 0.5 -0.5 0 0.5 -0.5 0 0.5 -0.5 0 0.5 -0.5 0 0.5 -0.5 0 0.5 -0.5 0 0.5 -0.5 0 0.5 -0.5 0 0.5 -0.5 0 0.5 -0.5 0 0.5 -0.5 0 0.5 -0.5 0 0.5 -0.5 0 0.5 |P | = 13 |P | = 21 |P | = 34 |P | = 55 |P | = 89 |P | = 144 |P | = 233 |P | = 377 Figure 5.3: For the set K (red), K0(P ) (green/light) under-approximates V iab[0,1](K,U) (outlined in thick black lines via [87]) using Al- gorithm 5.1 under the double integrator dynamics. A finer time interval partition results in better approximation. 13 21 34 55 89 144 233 377 100 150 200 250 300 350 400 450 E rr or (# of gr id p oi n ts ) |P | 8.7% 6.1% 4.2% 3.6% 3.0% 3.0% 2.6% 2.4% Figure 5.4: Convergence plot of the error as a function of |P | for the double- integrator example. Error is quantified as the fraction of grid points (total of 71×71) contained in the set difference between the level set approximation of the viability kernel and its piece- wise ellipsoidal under-approximation. 95 5.2. Computational Algorithms this under-approximation may be too conservative—particularly when the constraints K and U are nonsymmetric with respect to the origin. In such a case we would typically have to reevaluate our time interval partitioning and choose a P with extremely short sub-intervals so that it offsets the potential conservatism. This is clearly not an ideal approach. Instead consider the respective linear differential inclusion ẋ(t) ∈ AK ⊕BU . (5.22) Notice that the set AK ⊕ BU is a compact (closed and bounded) subset of X since both K and U are compact. With this in mind, a possibly less conservative approach to construct K↓(P ) may be as follows. 1. Over-approximate the set AK ⊕ BU by computing tight enclosing el- lipsoids over a finite subset of directionsM⊂ {l ∈ Rn | 〈l, l〉 = 1} and choose the one with the smallest volume: E(σ,Σ) := minvol l∈M ( AK E+l⊕ BU ) . (5.23) Here we have used E+l⊕ to denote the tight enclosing ellipsoidal approx- imation of the Minkowski sum in the direction l. 2. Multiply this set by a factor of ‖P‖: ‖P‖E(σ,Σ) = E(‖P‖σ, ‖P‖2Σ). (5.24) 3. Under-approximate the set K E(‖P‖σ, ‖P‖2Σ) by computing tight enclosed ellipsoids over directions in M and choose the one with the largest volume to obtain K↓(P ): K↓(P ) = maxvol l∈M ( K E−l E(‖P‖σ, ‖P‖2Σ) ) . (5.25) Here the operator E−l returns the tight enclosed ellipsoidal approxi- 96 5.3. Practical Examples mation of the Pontryagin difference (or the erosion operator) between its operands in the direction l. All of the above steps can be easily and efficiently performed in ET. 5.3 Practical Examples All computations are performed on a dual core Intel-based computer with 2.8 GHz CPU, 6 MB of L2 cache and 3 GB of RAM running single-threaded 32-bit Matlab 7.5. 5.3.1 Flight Envelope Protection (Continuous-Time) Consider the linearized longitudinal aircraft dynamics ẋ(t) = Ax(t)+Bδe(t), A = −0.003 0.039 0 −0.322 −0.065 −0.319 7.740 0 0.020 −0.101 −0.429 0 0 0 1 0 , B = 0.010 −0.180 −1.160 0 with state x = [u, v, θ̇, θ]T ∈ R4 comprised of deviations in aircraft velocity [ft/s] along and perpendicular to body axis, pitch-rate [crad/s],3 and pitch angle [crad] respectively, and with input δe ∈ [−13.3◦, 13.3◦] ⊆ R the elevator deflection. These matrices represent stability derivatives of a Boeing 747 cruising at an altitude of 40 kft with speed 774 ft/s [16]. The state constraint set K = E 0 0 2.18 0 , 1075.84 0 0 0 0 67.24 0 0 0 0 42.7716 0 0 0 0 76.0384 represents an ellipsoidal (under-)approximation of the flight envelope. We require x(t) ∈ K ∀t ∈ [0, 2]. 3crad = 0.01 rad ≈ 0.57◦. 97 5.3. Practical Examples Figure 5.5: 3D projections of the under-approximation of V iab[0,2](K,U) for Example 5.3.1. The flight envelope K is the red/light trans- parent region. The green/dark piecewise ellipsoidal sets under- approximate the viability kernel. A uniform partition P is chosen such that |P | = 400. Algorithm 5.1 (with |M| = 8 consisting of the standard basis vectors in R4 and four additional randomly generated vectors) computes via ET a piecewise ellipsoidal under- approximation of the viability kernel V iab[0,2](K,U) as shown in Figures 5.5 and 5.6. Note that for any state belonging to this set, there exists an input that can protect the flight envelope over the specified time horizon. The overall computation time was roughly 10 mins. In comparison, the level set approximation of the viability kernel (also shown in Figure 5.6) is computed in 5.4 hrs with significantly larger memory footprint over a grid with 45 nodes in each dimension (still a rather coarse grid) using the Level Set Toolbox [87]. Since the computed sets are 4D, we plot a series of 3D and 2D projections of these 4D objects. 98 5.3. Practical Examples -20 0 20 -10 -5 0 5 10 -20 0 20 -10 -5 0 5 10 -20 0 20 -10 -5 0 5 10 -10 0 10 -10 -5 0 5 10 -10 0 10 -10 -5 0 5 10 -10 0 10 -10 -5 0 5 10 x1 x1 x1 x3x2x2 x 2 x 3 x 3 x 4 x 4 x 4 Figure 5.6: 2D projections of the under-approximation of V iab[0,2](K,U) for Example 5.3.1. The constraint set K (red/dark) and a piecewise ellipsoidal under-approximation of the viability ker- nel (green/light) are shown. The level set approximation of the viability kernel, computed via [87], is outlined in thick black lines. 99 5.3. Practical Examples 5.3.2 Safety in Anesthesia Automation (Discrete-Time) To improve patient recovery, lessen anesthetic drug usage, and reduce time spent at drug saturation levels, a variety of approaches to controlling depth of anesthesia have been proposed e.g. in [10, 26, 52, 85, 97, 109, 115]. Over the past few years, an interdisciplinary team of researchers at the University of British Columbia has been developing an automated drug de- livery system for anesthesia. As part of this effort, an open-loop bolus-based neuromuscular blockade system was developed and clinically validated in [34]. Discrete-time Laguerre-based LTI models of the dynamic response to rocuronium were identified using data collected from more than 80 patients via clinical trials.4 To obtain regulatory certificates to fully close the loop while employing an infusion-based administration of the drugs, mathemati- cal guarantees of safety of the system are likely to be required. The viability kernel can provide such guarantees. Let us consider the problem of comput- ing the viability kernel for a constrained dynamical system that represents the pharmacological response of a patient under anesthesia subject to ther- apeutic bounds. Patient Model and Constraints Consider the following discrete-time LTI system describing the Laguerre dynamics of a patient [34, 35]: x(t+ 1) = Ax(t) +Bu(t), (5.26) y(t) = Cx(t) (5.27) with time step t ∈ Z+, state vector x(t) ∈ R6, input (rocuronium infusion rate [mg/kg/min]) u(t) ∈ R, and output (pseudo-occupancy, a metric related to the patient’s plasma concentration of anesthetic e.g. rocuronium) y(t) ∈ 4The states of such a model correspond to Laguerre polynomials whose weighted sum attempts to closely rebuild the shape of the actual impulse (bolus) response of the patient; See [122] for more details on the Laguerre modeling framework. 100 5.3. Practical Examples R. The sampling interval is 20 s and the system matrices are: A = 0.9960 0 0 0 0 0 0.0080 0.9960 0 0 0 0 −0.0080 0.0080 0.9960 0 0 0 0.0079 −0.0080 0.0080 0.9960 0 0 −0.0079 0.0079 −0.0080 0.0080 0.9960 0 0.0079 −0.0079 0.0079 −0.0080 0.0080 0.9960 , B = [ 0.0894 −0.0890 0.0886 −0.0883 0.0879 −0.0876 ]T , C = [ 18.5000 8.2300 3.5300 4.3400 3.7000 3.0700 ] . The constraint set (desired clinical effect) is specified in terms of the pseudo- occupancy level and the input is bounded above and below by hard physical constraints: y(t) ∈ K0 := [0.1, 1],u(t) ∈ U0 := [0, 0.8]. (5.28) Reformulating the Problem Note that the therapeutic constraint is specified in the output space (as opposed to the state space) and the output signal y should track a reference setpoint that lies within K0. To perform our desired analysis on this system, we reformulate the problem by projecting the output bounds onto the state space while making the control action regulatory. For brevity we drop the time argument from the state, input, and output notations. Projection of Bounds onto the State Space Consider the (nonsingu- lar) linear transformation[ C 05×1 I5 ] [ x1 x2 · · · x6 ]T =: [ w1 w2 · · · w6 ]T . The states w2, . . . , w6 are the Laguerre states x2, . . . , x6. In the new co- ordinate space, the bounds are state space constraints on the first state 101 5.3. Practical Examples w1 := Cx = y. Tracking vs. Regulating We perform an affine change of coordinates and shift the equilibrium point to the origin. This is done by augmenting the state vector in the w-space with the reference output signal y∗ and applying a basis translation so that in the new coordinates the first state w1 := y becomes z1 := y − y∗: 1 0 · · · 0 −1 0 1 . . . 0 ... . . . . . . . . . ... ... . . . 1 0 0 · · · · · · 0 1 y w2 ... w6 y∗ = y − y∗ w2 ... w6 y∗ =: z1 z2 ... z6 z7 . (5.29) Let u(t) = uss be the steady state control input needed for tracking a con- stant setpoint y∗(t) = y∗ = 0.9. This value can be easily calculated using a standard state-feedback procedure from[ xss uss ] = [ A− I B C 0 ]−1 [ 0 y∗ ] , (5.30) where xss denotes the steady state equilibrium. To complete the reformula- tion, we deduct uss from the control set of the transformed system. The new constraints for the transformed, extended system z(t + 1) = Ãz(t) + B̃u(t), y(t) = C̃z(t), with Ã ∈ R7×7, B̃ ∈ R7×1, C̃ ∈ R1×7, are as follows: z(t) ∈ K := (K0 − y∗)× R6,u(t) ∈ U := U0 − uss. (5.31) In this new seven-dimensional state space the first state z1 represents the drug pseudo-occupancy minus its setpoint value of 0.9 units, the next five states are the second to sixth Laguerre states transformed from the original coordinates, and the last state z7 is a constant corresponding to the pseudo-occupancy setpoint. The states are assumed to be constrained 102 5.4. Summary and Further Discussions by a slab in R7 that is only bounded in the z1 direction. Note that with this formulation, the last state z7 is allowed to take on values that are not needed; of actual interest is the behavior of the remaining states when z7 equals the pseudo-occupancy setpoint y∗. The input constraint set is a one-dimensional ellipsoid. To under-approx- imate the state constraint with a non-degenerate ellipsoid we use a priori knowledge about the typical values of the (Laguerre) states z2, . . . , z6 and bound them by an ellipsoid with a large spectral radius of λmax = 30 in those directions. (This imposed constraint can be further relaxed if necessary.) Guaranteeing that this ellipsoidal target set K, which is our desired clinical effect, is not violated during the surgery provides a certificate of safety of the closed-loop system. Therefore, for a 30 min surgery for instance, we require z(t) ∈ K ∀t ∈ [0, 90] despite bounded input authority. Using appropriately synthesized infusion policies, the states belonging to the viability kernel of K under the extended system will never leave the desired clinical effect for the duration of the surgery. We under-approximate V iab[0,90](K,U) in 986 s using Algorithm 5.2 with |M| = 30. Of the 30 randomly chosen initial directions used in the ellipsoidal computations, 15 resulted in nonempty ellipsoids that make up the piecewise ellipsoidal under-approximation of the viability kernel (Figure 5.7). Note that no similar computations are currently possible in such high dimensions using Eulerian methods directly. 5.4 Summary and Further Discussions In this chapter we presented a new connection between the viability kernel (and by duality, the minimal reachable tube) and the maximal reachable sets of possibly nonlinear systems. Owing to this connection, the efficient and scalable Lagrangian techniques that were traditionally developed for the approximation of the maximal reachability constructs can now be used to ap- proximate the viability kernel. Motivated by a high-dimensional problem of guaranteed safety in control of anesthesia, we proposed a scalable algorithm that computes a piecewise ellipsoidal under-approximation of the viability 103 5.4. Summary and Further Discussions -1 0 1 -20 0 20 -1 0 1 -20 0 20 -1 0 1 -20 0 20 -1 0 1 -20 0 20 -1 0 1 -20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20 -20 0 20 z1 z1 z1 z1 z1 z2 z2 z2 z2 z3 z3 z3 z4 z4 z5 z 2 z 3 z 4 z 5 z 6 z 3 z 4 z 5 z 6 z 4 z 5 z 6 z 5 z 6 z 6 Figure 5.7: 2D projections of the under-approximation of V iab[0,90](K,U) for Example 5.3.2 for the first six states when z7 equals the setpoint value. The constraint set K (blue/dark) and a piece- wise ellipsoidal under-approximation of the provably safe regions (green/light) are shown. 104 5.4. Summary and Further Discussions kernel for LTI systems based on ellipsoidal techniques for reachability. Empirically quantifying the computational complexity of the piecewise ellipsoidal algorithm is a work under way for which we expect a polynomial complexity in the order of |M||P | (O(Rδ) +O(S)) where O(Rδ) is the com- plexity of computing the maximal reachable set along a given direction over the time interval δ and O(S) is the complexity of solving the SDP (5.19). Our preliminary tests for a chain of double integrators using a single direc- tion (yielding one ellipsoid in Rn) show that the technique scales relatively well up to about 35 dimensions, which is computed in less than 200 s. The presented connection between the viability kernel and the maximal reachable sets paves the way to synthesizing safety-preserving optimal con- trol laws in a more efficient and scalable manner. This scalable synthesis methodology as well as an extension of the results of this chapter to the dif- ferential games setting and the approximation of the discriminating kernel are presented in the following chapter. 105 Chapter 6 Robust Approximation and Scalable Safety-Preserving Control Synthesis In this chapter we will first extend the results presented in the previous chapter to the differential games setting—the case in which the system is perturbed by an unknown but bounded disturbance input or uncertainty in the model. We will then propose a safety-preserving control strategy based on the piecewise ellipsoidal algorithm discussed in Section 5.2. Although the synthesis of safety-preserving controllers is by no means a new research di- rection (cf. e.g. [5, 11, 81, 116]), the results presented here provide a scalable synthesis of such controllers for a class of constrained systems. We will make use of the definitions and preliminaries in Sections 2.5 and 2.6. We focus only on the continuous-time case. In the discrete-time case the Isaac’s condition does not in general hold and consequently the discussion is more involved. 6.1 The Discriminating Kernel in Terms of Robust Maximal Reachable Sets For the constrained system (2.29), i.e. ẋ(t) = f(x(t), u(t), v(t)), x(0) = x0 (6.1) 106 6.1. Discriminating Kernel vs. Robust Maximal Reach Sets with control input u(t) ∈ U and disturbance input v(t) ∈ V, we show that we can approximate Disc[0,τ ](K,U ,V) by considering a nested sequence of sets that are robustly reachable in small sub-intervals of [0, τ ]. Similar to before, we say that the vector field f is bounded on K if there exists a norm ‖·‖ : X → R+ and a real number M > 0 such that for all x ∈ K, u ∈ U , and v ∈ V we have ‖f(x, u, v)‖ ≤ M . If K is compact, every continuous vector field f is bounded on K. 6.1.1 Under-Approximating the Discriminating Kernel The discriminating kernel can be under-approximated via the recursion K|P |(P ) = K↓(P ), (6.2a) Kk−1(P ) = K↓(P ) ∩Reach]tk−tk−1(Kk(P ),U ,V) for k ∈ {1, . . . , |P |} (6.2b) where K↓(P ) := {x ∈ K | dist(x,Kc) ≥M‖P‖} is a subset of K deliberately chosen at a distance of M‖P‖ from its boundary. Proposition 6.1. For any partition P ∈ P([0, τ ]) the final set K0(P ) de- fined by the recurrence relation (6.2) satisfies K0(P ) ⊆ Disc[0,τ ](K,U ,V). (6.3) Proof. The proof is a straightforward extension of the proof of Proposi- tion 5.1: Since f is bounded on K, there exists a norm ‖·‖ and a real num- ber M > 0 such that ‖f(x, u, v)‖ ≤ M ∀(x, u, v) ∈ K × U × V. Now, fix a partition P of [0, τ ] and take a point x0 ∈ K0(P ). By the con- struction of K0(P ), this means that for each k = 1, . . . , |P | there is some point xk ∈ Kk(P ) and a control uk : [0, tk − tk−1] → U for any (non- anticipative) disturbance vk = ρ[uk] : [0, tk − tk−1] → V such that xk can be reached from xk−1 at time tk − tk−1 using uk. Thus, taking the con- catenation of the inputs uk and vk, we get a control input u : [0, τ ] → U for every disturbance input v = ρ[u] : [0, τ ] → V such that the solution 107 6.1. Discriminating Kernel vs. Robust Maximal Reach Sets x : [0, τ ]→ X to the initial value problem ẋ = f(x, u, v), x(0) = x0, satisfies x(tk) = xk ∈ Kk(P ) ⊆ {x ∈ K | dist(x,Kc) ≥ M‖P‖}. We claim that this guarantees that x(t) ∈ K ∀t ∈ [0, τ ]. Indeed, any t ∈ [0, τ) lies in some interval [tk, tk+1). Since f is bounded by M , we have ‖x(t)− x(tk)‖ ≤ ∫ t tk ‖ẋ(s)‖ds ≤M(t− tk) < M(tk+1 − tk) ≤M‖P‖. (6.4) Further, x(tk) ∈ Kk(P ) implies dist(x(tk),Kc) ≥ M‖P‖. Combining these, we see that dist(x(t),Kc) ≥ dist(x(tk),Kc)− ‖x(t)− x(tk)‖ > M‖P‖ −M‖P‖ = 0 (6.5) and hence x(t) ∈ K. Thus, x0 ∈ Disc[0,τ ](K,U ,V). Corollary 6.1 (Intermediate Discriminating Kernels). For any partition P ∈ P([0, τ ]) and every k ∈ {1, . . . , |P |} the sets Kk−1(P ) defined by the recurrence relation (6.2) satisfy Kk−1(P ) ⊆ Disc[0,τ−tk−1](K,U ,V) (6.6) with K|P |(P ) ⊆ Disc[0,0](K,U ,V) = K. Similar to the viability kernel case, the approximation here can be made to be arbitrarily precise by choosing a sufficiently fine partition. Proposition 6.2. Suppose that the vector field f : X×U×V → X is bounded on a set K ⊆ X . Then we have Disc[0,τ ]( ◦K,U ,V) ⊆ ⋃ P∈P([0,τ ]) K0(P ) ⊆ Disc[0,τ ](K,U ,V). (6.7) 108 6.2. Computational Algorithm & Control Synthesis In particular, when K is open,⋃ P∈P([0,τ ]) K0(P ) = Disc[0,τ ](K,U ,V). (6.8) Proof. The proof is a straightforward extension of the proof of Proposi- tion 5.2: The second inclusion in (6.7) follows directly from Proposition 6.1. To prove the first inclusion, take a state x0 ∈ Disc[0,τ ]( ◦K,U ,V). There exists a control u : [0, τ ] → U for every disturbance v = ρ[u] : [0, τ ] → V such that the trajectory x(·) satisfies x(t) ∈ ◦K ∀t ∈ [0, τ ]. Since ◦K is open, ∀x ∈ ◦K dist(x,Kc) > 0. Further, x : [0, τ ] → X is continuous so the func- tion t 7→ dist(x(t),Kc) is continuous on the compact set [0, τ ]. Thus, we can define d > 0 to be its minimum value. Take a partition P of [0, τ ] such that ‖P‖ < d/M . We need to show that x0 ∈ K0(P ). First note that our partition P is chosen such that dist(x(t),Kc) > M‖P‖ ∀t ∈ [0, τ ]. Hence x(tk) ∈ K|P |(P ) for all k = 0, . . . , |P |. To show that x(tk−1) ∈ Reach]tk−tk−1(Kk(P ),U ,V) for all k = 1, . . . , |P |, consider the to- kenizations {uk}k and {vk}k of the control u and the disturbance v, respec- tively, corresponding to P . It is easy to verify that for all k we can reach x(tk) from x(tk−1) at time tk−tk−1 using the control input uk for every vk = ρ[uk]. Thus, in particular, we have x0 = x(t0) ∈ Reach]t1−t0(K1(P ),U ,V). So x0 ∈ K0(P ). Hence Disc[0,τ ]( ◦K,U ,V) ⊆ ⋃P∈P([0,τ ])K0(P ). By approximatingDisc[0,τ ](K,U ,V) using the sub-interval maximal reach- able sets via recursion (6.2) we enable the use of scalable Lagrangian tech- niques for the computation of the discriminating kernel for high-dimensional systems. As we shall see, this will also allow us to compute the safety- preserving control laws more efficiently. 6.2 Computational Algorithm & Safety-Preserving Control Strategy The following corollary forms the basis of our safety-preserving feedback strategy. It follows directly from Proposition 6.1, its proof, and the recur- 109 6.2. Computational Algorithm & Control Synthesis sion (6.2). Corollary 6.2. For a fixed time partition P ∈P([0, τ ]) suppose K0(P ) 6= ∅. For any initial condition x0 ∈ K0(P ) the concatenation u of sub-interval control inputs uk : [0, tk− tk−1]→ U corresponding to robust maximal reach- able sets Reach]tk−tk−1(Kk(P ),U ,V) for k = 1, . . . , |P | is a safety-preserving control law that keeps the trajectory x(·) of the system (6.1) with x(0) = x0 contained in K for every v(t) = ρ[u](t) ∈ V for all time t ∈ [0, τ ]. Notice that via Corollary 6.2 any Lagrangian method that facilitates the synthesis of maximal reachability controllers can be employed to form a safety-preserving policy. One such Lagrangian method is the ellipsoidal techniques [70] implemented in Ellipsoidal Toolbox (ET) [73]. Consider the case in which (2.1) is a linear time-invariant system ẋ(t) = Ax(t) +Bu(t) +Gv(t) (6.9) with A ∈ Rn×n, B ∈ Rn×mu , and G ∈ Rn×mv . We will further assume that the constraints K, U , and V are (or can be closely under-approximated by) nonempty compact ellipsoids. In Chapter 5 we introduced a scalable piecewise ellipsoidal algorithm (Al- gorithm 5.1) based on the ellipsoidal techniques for under-approximating the viability kernel V iab[0,τ ](K,U). That result can easily be extended so that the generated set approximates the discriminating kernel Disc[0,τ ](K,U ,V) instead by simply computing the intermediate maximal reachable sets for adversarial inputs in line with the analysis in Proposition 6.1. Before we describe the corresponding safety-preserving feedback strategy, let us sum- marize the aforementioned algorithm used to under-approximate the dis- criminating kernel. 6.2.1 Piecewise Ellipsoidal Approximation of the Discriminating Kernel Similar to the discussions of Section 5.2.1, the robust maximal reachable set computed using the ellipsoidal techniques for a single direction is an 110 6.2. Computational Algorithm & Control Synthesis ellipsoidal subset of the actual robust maximal reachable set. With this in mind, and recalling that the function #(·) maps a set to its maximum volume inscribed ellipsoid, we can state the following. Proposition 6.3. For a fixed partition P ∈P([0, τ ]) and a terminal direc- tion `τ ∈M, the recursion K ∗[`τ ] k−1 = #(K|P |(P ) ∩Reach][`τ ]tk−tk−1(K∗[`τ ]k (P ),U ,V)) for k ∈ {1, . . . , |P |} (6.10) with K ∗[`τ ] |P | (P ) = K|P |(P ) defined as in (6.2a) generates an ellipsoidal set K ∗[`τ ] 0 (P ) such that⋃ `τ∈M K ∗[`τ ] 0 (P ) := K ∗ 0 (P ) ⊆ Disc[0,τ ](K,U ,V). (6.11) Proof. The proof is similar to that of Proposition 5.4 which discusses the case in which V = {0}. The set K∗0 (P ) is therefore a piecewise ellipsoidal under-approximation of the discriminating kernel. Notice however that the final #(·) operation when k = 1 is not necessary if a closed-form piecewise ellipsoidal expression is not needed. Indeed, one can easily verify that K∗0 (P ) ⊆ K|P |(P ) ∩ ⋃ `τ∈M Reach ][`τ ] t1 (K ∗[`τ ] 1 (P ),U ,V) ⊆ Disc[0,τ ](K,U ,V). (6.12) Finally, note that similar to (6.6) the intermediate discriminating kernels are under-approximated via⋃ `τ∈M K ∗[`τ ] k−1 (P ) := K ∗ k−1(P ) ⊆ Disc[0,τ−tk−1](K,U ,V). (6.13) Remark 6.1. The under-approximation K|P |(P ) = K↓(P ) may be con- structed by either eroding the set K by a ball of radius M‖P‖ (for a given 111 6.2. Computational Algorithm & Control Synthesis uniform bound M), or alternatively, by using the method described on page 94 with the difference that the right-hand side of (5.23) is now replaced by minvol l∈M ( E+l⊕ {AK, BU , GV} ) . (6.14) 6.2.2 Safety-Preserving Feedback Policy Suppose U := E(µ,U) (6.15) with center µ ∈ Rmu and shape matrix U ∈ Rmu×mu . Based on the piecewise ellipsoidal algorithm described above, we can form a control policy taking values pointwise on U that keeps the trajectory of the system in K over the entire time horizon (despite the actions of the disturbance) by concatenat- ing the sub-interval robust maximal reachability control laws according to Corollary 6.2. Suppose that in computing the reachable setReach ][`τ ] tk−tk−1(K ∗[`τ ] k (P ),U ,V) we can also compute Reach ][`τ ] tk−t(K ∗[`τ ] k (P ),U ,V) ∀t ∈ (tk−1, tk). (6.16) That is, we compute the entire reachable tube over the kth sub-interval. For fixed `τ ∈M and k ∈ {1, . . . , |P |} let xc[`τ ]k (t−tk−1) and X−[`τ ]`,k (t−tk−1) de- note the center and the shape matrix of the ellipsoidReach ][`τ ] tk−t(K ∗[`τ ] k (P ),U ,V) at time t ∈ [tk−1, tk] with K∗[`τ ]k (P ) = E(xc[`τ ]k (tk − tk−1), X−[`τ ]`,k (tk − tk−1)) and the terminal direction `τ . Define a shorthand notation R(t, k) := ⋃ `τ∈M R[`τ ](t, k) := ⋃ `τ∈M Reach ][`τ ] tk−t(K ∗[`τ ] k (P ),U ,V) = ⋃ `τ∈M E(xc[`τ ]k (t− tk−1), X−[`τ ]`,k (t− tk−1)), (6.17) and suppose R(0, 1) 6= ∅. 112 6.2. Computational Algorithm & Control Synthesis qperf qsafe û(x(t), t) ∈ U û(x(t), t) x(t) ∈ ◦R[γ](t, k) x(t) 6∈ ◦R[γ](t, k) x(t) 6∈ ◦R(t, k) x(t) ∈ ◦R[γ](t, k) x(t) ∈ ◦R[¯̀τ ](t, k) ∃ ¯̀τ 6= γ s.t.γ := ¯̀τ γ := ¯̀τ x(t) ∈ ◦R[¯̀τ ](t, k) ∃ ¯̀τ 6= γ s.t. γ := `τ x0 ∈ R[`τ ](0, 1) ∃`τ ∈M s.t. = ψBU(l[γ]) Figure 6.1: The graph of the hybrid automaton H representing the safety- preserving controller. Now, consider a controller described by the hybrid automaton H = (Q,Σe,Σi, Init,Dom,E,G,R,Ufb) (6.18) where Q = {qperf, qsafe} is the set of discrete states with qperf representing the case in which the controller is free to choose any value in U (“per- formance” mode) and with qsafe representing the case in which the con- troller is required to return the optimal safety-preserving law (“safety” mode); See Figure 6.1. The inputs to the controller are drawn from the sets Σe ⊆ X × [0, τ ] (external input) and Σi ⊆ M (internal input). The external input is the pair (x(t), t) ∈ Σe and the internal input is the direc- tion vector γ ∈ Σi. The initial state of the automaton Init ⊆ Q× Σe × Σi is assumed to be Init = (qperf, x0, 0, {γ | ∃`τ ∈ M, x0 ∈ R[`τ ](0, 1), γ = `τ}). The domains Dom(·, γ, t) : Q → 2X of the automaton for every γ ∈ Σi and t ∈ [0, τ ] are Dom(qperf, γ, t) = ◦R[γ](t, k) and Dom(qsafe, γ, t) = (Dom(qperf, γ, t)) c. The domains specify (γ, t)-varying invariants for every (x(t), t) ∈ Σe that must be satisfied in each mode. The edges E ⊆ Q×Q are 113 6.2. Computational Algorithm & Control Synthesis E = {(qperf, qperf), (qperf, qsafe), (qsafe, qperf)1, (qsafe, qperf)2} where subscripts are used when necessary in order to distinguish between the two edges of the automaton that enable transitions from qsafe to qperf that possess dif- ferent properties (e.g. guards and reset rules). The guards G(·, γ, t) : E → 2X for every γ ∈ Σi and t ∈ [0, τ ] are conditions on (x(t), t) ∈ Σe de- fined as: G(qperf, qsafe, γ, t) = ( ◦R(t, k))c, G((qsafe, qperf)1, γ, t) = ◦R[γ](t, k), and G((qsafe, qperf)2, γ, t) = G(qperf, qperf, γ, t) = { ◦R[¯̀τ ](t, k) for some ¯̀τ ∈ Σi, ¯̀τ 6= γ}. The domains and the guards are chosen in terms of the in- terior of the set R[γ](t, k) to ensure that the automaton H is non-blocking and that transitions over E can take place when necessary. A transition corresponding to an edge is enabled for every t if x(t) satisfies its guard. For example, the automaton can make a transition from qsafe to qperf over the edge (qsafe, qperf)2 ∈ E for a fixed (γ, x(t), t) ∈ Σi × Σe if ∃¯̀τ 6= γ s.t. x(t) ∈ R[¯̀τ ](t, k). The map R : E × Σi → Σi resets the internal input via R((qsafe, qperf)2, γ) = ¯̀τ , R(qperf, qperf, γ) = ¯̀τ , and R((qsafe, qperf)1, γ) = R(qperf, qsafe, γ) = γ. Finally, the output of H is a set-valued map Ufb : Q× Σi × Σe → 2U given byUfb(qperf, γ, x(t), t) = U ;Ufb(qsafe, γ, x(t), t) = ψBU (l[γ]), (6.19) where ψBU : Rn → U , l[γ] 7→ µ− UBTl[γ]〈l[γ], BUBTl[γ]〉−1/2 (6.20) is chosen so thatBψBU (l[γ]) is the support vector of the setBU = E(Bµ,BUBT) ⊆ X in the direction −l[γ] ∈ Rn with l[γ] = l[γ](x(t), t) = (X −[γ] `,k (t− tk−1))−1(x(t)− xc[γ]k (t− tk−1)). (6.21) (This strategy is based on the optimal control design presented in [69] and 114 6.2. Computational Algorithm & Control Synthesis [75].) Notice that we allow non-determinism in the mode transitions of the hybrid automaton to formulate a non-restrictive policy (in the sense of [116]); q = qsafe only when safety is at stake, and q = qperf otherwise. As we shall see, the primary objective of the controller, i.e. preserving safety, is achievable regardless of this behavior. Theorem 6.1 (Safety-Preserving Controller). For a given partition P ∈ P([0, τ ]) for any x0 ∈ K∗0 (P ) where K∗0 (P ) is the piecewise ellipsoidal set obtained through (6.10)–(6.11), the feedback policy u(t) = û(x(t), t) ∈ Ufb (6.22) generated by the hybrid automaton H keeps the trajectory x of the system (6.9) with initial condition x(0) = x0 contained in K for any disturbance v(t) = ρ[u](t) ∈ V for all time t ∈ [0, τ ]. Proof (Adapted from [69]). Let k ∈ {1, . . . , |P |} be the unique integer such that t ∈ [tk−1, tk) =: Tk. We prove that safety is preserved in each mode for any given (γ, x(t), t) ∈M×X × Tk for all v(t) ∈ V. First, note that if for every k, x(tk−1) ∈ K∗k−1(P ) (which, as we shall see, is indeed the case) then x(tk−1) ∈ R(tk−1, k) since K∗k−1(P ) ⊆ R(tk−1, k). Fix γ = `τ ∈M and k and let x(tk−1) ∈ R[γ](tk−1, k). Define a continuously differentiable value function V [γ] k : X × Tk → R such that V [γ] k (x(t), t) = dist 2(x(t),R[γ](t, k)). (6.23) Notice that V [γ] k (x(t), t) = 0 for x(t) ∈ ◦R[γ](t, k) and V [γ]k (x(t), t) ≥ 0 for x(t) 6∈ ◦R[γ](t, k). We use the convention R[γ](tk−1, k) ≡ { x ∈ X | V [γ]k (x, tk−1) ≤ 0 } . (6.24) and assume without loss of generality that the disturbance plays a worst-case game against the control. Consider the case in which x(t) 6∈ ◦R[γ](t, k) and the active mode of 115 6.2. Computational Algorithm & Control Synthesis the automaton is qsafe. Clearly, we have that dist(x(t),R[γ](t, k)) ≥ 0. As shown in [69], computing the Lie derivative of V [γ] k (x, t) along the system (6.9) yields d dt V [γ] k (x(t), t) = 2 dist(x(t),R[γ](t, k)) × d dt dist(x(t),R[γ](t, k)) (6.25) with d dt dist(x(t),R[γ](t, k)) = ∂ ∂t dist(x(t),R[γ](t, k)) + 〈 ∂ ∂x dist(x(t),R[γ](t, k)), ẋ(t) 〉 (6.26) = 〈 l[γ], Bu(t) +Gv(t) 〉 + ρBU (−l[γ]) − 〈l[γ], Gv(t)〉 (6.27) where ρBU (l) := supz∈BU 〈l, z〉 is the support function of the set BU in the direction l ∈ Rn and l[γ] is the unique direction vector given by (6.21) [75]. We refer the reader to [69, Sections 1.4 and 1.8] for additional details on how (6.27) is obtained from (6.26). Notice that the right-hand side simplifies to〈 l[γ], Bu(t) 〉 + ρBU (−l[γ]) = 〈 l[γ], Bu(t) 〉 + sup z∈BU 〈− l[γ], z〉 = 〈 l[γ], Bu(t) 〉− inf z∈BU 〈 l[γ], z 〉 . (6.28) Recalling the fact that the set BU is a compact ellipsoid and invoking (6.20), one can verify that d dt dist(x(t),R[γ](t, k)) = 0 for u(t) = ψBU (l[γ]); (6.29) d dt dist(x(t),R[γ](t, k)) ≥ 0 for any u(t) ∈ U . (6.30) When the disturbance does not play a worst-case game, (6.27) is replaced by ddt dist(x(t),R[γ](t, k)) ≤ 〈l[γ], Bu(t) + Gv(t)〉 + ρBU (−l[γ]) − ρGV(l[γ]) 116 6.2. Computational Algorithm & Control Synthesis which would imply ddt dist(x(t),R[γ](t, k)) ≤ 0 for u(t) = ψBU (l[γ]) and ∀v(t) = ρ[u](t) ∈ V. On the other hand, when x(t) ∈ ◦R[γ](t, k) and the active mode of the controller is qperf we have that dist(x(t),R[γ](t, k)) = 0. Thus the derivative of the distance function need not be examined as the Lie derivative of the value function (6.25) is automatically zero regardless of the value of u(t). Combining the above we see that with u(t) = û(x(t), t) ∈ Ufb(q, γ, x(t), t) we have d dt V [γ] k (x, t) ∣∣∣ (6.9) ≤ 0 ∀(q, γ, x, t, v) ∈ Q×M×X × Tk × V. (6.31) This yields∫ t tk−1 d ds V [γ] k (x(s), s)ds = V [γ] k (x(t), t)− V [γ]k (x(tk−1), tk−1) ≤ 0 (6.32) which in turn implies (via (6.24)) V [γ] k (x(t), t) ≤ V [γ]k (x(tk−1), tk−1) ≤ 0 ∀t ∈ Tk. (6.33) Therefore, for every solution x of the differential inclusion ẋ(t) ∈ Ax(t) ⊕ BUfb(q, γ, x(t), t)⊕GV we will have x(t) ∈ R[γ](t, k) = Reach][γ]tk−t(K ∗[γ] k (P ),U ,V) ∀t ∈ Tk, (6.34) x(tk) ∈ K∗[γ]k (P ). (6.35) Notice that according to (6.13) K ∗[γ] k (P ) ⊆ ⋃ γ∈M K ∗[γ] k (P ) := K ∗ k(P ) ⊆ Disc[0,τ−tk](K,U ,V) ⊆ K. (6.36) Therefore, x(tk) ∈ K, for every k ∈ {1, . . . , |P |} granted that the feedback control law (6.22) is applied. 117 6.2. Computational Algorithm & Control Synthesis It remains to show that x(t) ∈ K for all t ∈ (tk−1, tk). Indeed, by construction of the recursion (6.10) via Proposition 6.1 we have that for all t ∈ Tk and k ∈ {1, . . . , |P |} Reach ][γ] tk−t(K ∗[γ] k (P ),U ,V) ⊆ {x ∈ K | dist(x,Kc) ≥M‖P‖}. (6.37) Therefore combining all this we conclude that with u(t) ∈ Ufb and x0 ∈ R(0, 1) we have xu,ρ[u]x0 (t) ∈ K ∀t ∈ [0, τ ] for every v(t) = ρ[u](t) ∈ V. 6.2.3 Remarks and Practical Considerations Respecting the Intermediate Kernels The strategy (6.22) will ensure that the trajectories of the closed-loop system evolve in the interior or else on the boundary of R(t, k) for all t ∈ [0, τ ] and k ∈ {1, . . . , |P |}. Thus as we have seen in (6.35) and (6.36), for a given P ∈ P([0, τ ]) even if the disturbance always plays a worst-case game the trajectories satisfy x(tk) ∈ K∗k(P ) ⊆ Disc[0,τ−tk](K,U ,V) for all k. Shared Boundary Points In practice for common points that are on the boundaries of two or more el- lipsoids Reach ][`τ ] tk−t(K ∗[`τ ] k (P ),U ,V) at every time instant t, any one of these ellipsoids can be used for the computation of l[`τ ](x(t), t) in (6.21) and its corresponding optimal control law ψBU (l[`τ ](x(t), t)) in mode qsafe. Any such control will ensure that the trajectory remains within the corresponding ro- bust maximal reachable tube while being steered towards K ∗[`τ ] k (P ) (which, as we have seen, will ultimately ensure constraint satisfaction and safety). To allow a more permissive control policy, the automaton H could be modified by adding a self-loop in mode qsafe (i.e. an edge (qsafe, qsafe)) so that among all possible safe control laws corresponding to those ellipsoids that share the common boundary point, one that simultaneously better satisfies a given performance criterion is selected. 118 6.2. Computational Algorithm & Control Synthesis Conservatism of K∗0 (P ): Synthesis for x0 6∈ K∗0 (P ) When the initial condition x0 lies inside of the true discriminating ker- nel but outside of its piecewise ellipsoidal under-approximation, i.e. x0 ∈ Disc[0,τ ](K,U ,V)\K∗0 (P ) for a given P ∈P([0, τ ]), the control policy (6.22) will not in general guarantee that the constraint K is respected for all t ∈ [0, τ ]. However, with the following modification the feedback law will attempt to keep the trajectory within K: For every t ∈ [0, s], s ∈ (0, τ ], with s being the first instant such that x(s) ∈ R(s, k′) for some k′ ∈ {1, . . . , |P |}, for every k ≤ k′ the direction vector l[`τ ](x(t), t) in (6.21) is modified to be the direction that corresponds to any `τ ∈ M for which dist(x(t), Reach ][`τ ] tk−t(K ∗[`τ ] k (P ),U ,V)) = min. We will demonstrate this us- ing an example in Section 6.3. Finally, note that for x0 6∈ Disc[0,τ ](K,U ,V) the modified control law described above may still be able to keep x in K over [0, τ ], only if the disturbance does not play optimally against the control. (In formulation of Disc[0,τ ](K,U ,V) we always assume that the disturbance plays its worst- case game.) On the other hand, in the deterministic case (V = {0}), for x0 6∈ V iab[0,τ ](K,U) there does not exist a control law that can keep the trajectory within K over [0, τ ]. Therefore any safety-preserving strategy ultimately fails. Smooth Transition Between Safety and Performance When the mode qperf of the controller is active, the strategy Ufb(qperf, γ, x(t), t) = U for some (γ, x(t), t) ∈ Σi×Σe is returned and any performance-satisfying control input in U can be chosen and applied to the system; denote this input as uperf(t). On the other hand, when the mode qsafe is active, the con- troller returns Ufb(qsafe, γ, x(t), t) = ψBU (l[γ]) and only this specific safety- preserving control law must be applied to the system to ensure safety; we denote this input as usafe(t). Choosing uperf arbitrarily without considering the main objective of the closed-loop system (preserving safety) may result in excessive switching be- tween the two modes of the controller (whose main priority is to preserve 119 6.2. Computational Algorithm & Control Synthesis safety). Thus the resulting control policy could have a high switching fre- quency and the controller could end up spending a significant amount of time at the extremum points of the input constraint set. Such a policy, among other shortcomings, may be hard on actuators in a practical setting. An ideal control policy in mode qperf should be a combination of both uperf and usafe (even though technically the safety component usafe is not needed when the controller is in mode qperf). One solution is a simple convex combination of these two inputs. That is, for every (γ, x(t), t) ∈ Σi × Σe and a given domain ◦R[γ](t, k) = ◦E(xc[γ]k (t− tk−1), X−[γ]`,k (t− tk−1)) in qperf we shall choose an input such that u(t) = (1− βα[φ[γ]](x(t), t))uperf(t) + βα[φ[γ]](x(t), t)usafe(t), (6.38) with βα[φ [γ]](x(t), t) := 1 if φ[γ](x(t), t) ≥ 1; 1 1−α(φ [γ](x(t), t)− α) if α ≤ φ[γ](x(t), t) ≤ 1; 0 if φ[γ](x(t), t) ≤ α, (6.39) where α ∈ [0, 1) is a design parameter (Figure 6.2) and φ[γ] : Σe → R+, (6.40) (x(t), t) 7→ 〈(x(t)− xc[γ]k (t− tk−1)), (X −[γ] `,k (t− tk−1))−1(x(t)− xc[γ]k (t− tk−1)) 〉 . (6.41) Note that in qperf for fixed γ and t we have that dom(φ [γ](·, t)) = ◦R[γ](t, k). Therefore, range(φ[γ](·, t)) = [0, 1). This is true simply because the set R[γ](t, k) is an ellipsoid, and therefore by definition, the one sub-level sets of the function φ[γ](·, t) in qperf form the interior of R[γ](t, k). That is, ◦R[γ](t, k) := {x ∈ X | φ[γ](x, t) < 1}. Notice that φ[γ](x(t), t) determines how “deep” inside the domain of qperf the trajectory is at time t by evaluating the one sub-level sets of the ellipsoid 120 6.2. Computational Algorithm & Control Synthesis α 1 1 0 φ β Figure 6.2: βα as a function of φ for a given design parameter α ∈ [0, 1). R[γ](t, k) (the closure of the domain of qperf) at x(t). The closer the trajec- tory is to the boundary of R[γ](t, k) the greater the value of βα[φ[γ]](x(t), t) and therefore the more pronounced the safety component usafe(t) of the in- put will be. On the other hand, if the trajectory is deep inside the domain of qperf, φ [γ](x(t), t) tends to α, and therefore βα[φ [γ]](x(t), t) goes to zero. As a result, a greater emphasis is given to the performance component uperf(t) of the input. As such, the parameter α determines where within the domain of qperf the safety component usafe(t) should kick in. Clearly, larger values of α imply more emphasis on performance, while smaller values yield smoother transition between qperf and qsafe. Note that the limit of the control law u(t) in (6.38) as x(t)→ ∂R[γ](t, k) is usafe(t). Therefore the control is continuous across the automaton’s transition to qsafe. Such a policy will ensure a gradual change of the effective component of the control law from one form to the other, resulting in less switching frequency between performance and safety. We will show this policy using a number of examples in Section 6.3. Intermediate Maximal Reachable Sets vs. Tubes Notice that while in the approximation of the discriminating/viability ker- nel via recursion (6.10) only the final intermediate maximal reachable sets Reach ][`τ ] tk−tk−1(K ∗[`τ ] k (P ),U ,V) are used, in the control synthesis scheme de- scribed above (which is based on the ellipsoidal techniques [69]) to form a 121 6.3. Numerical Examples proper state-feedback law all intermediate reachable tubes⋃ t∈Tk Reach ][`τ ] tk−t(K ∗[`τ ] k (P ),U ,V) (6.42) are required and must be recorded. Chattering and Zeno Executions In theory, the continuous-time evolution of the dynamics in each mode of H and the particular formulation of the automaton may cause chattering (fre- quent switches between the two modes) or Zeno behavior (infinite switches in finite time). While such undesirable phenomena may be avoided by impos- ing certain conditions such as dwell-time [92] (requiring that the automaton remains in each mode for a non-zero amount of time) or other techniques as discussed in e.g. [2, 29, 53], in practice for a machine implementation of our control algorithm a finite number of sets can be used to form the intermediate maximal reachable tubes due to the discrete nature of numer- ical evaluations of the integration of the differential equations. This, in principle, imposes a constant dwell-time (which equals to the integration time-step) and thus ensures that the automaton is well-behaved. In addi- tion, the previously presented scheme for smooth transition between safety and performance (page 119) can also prevent chattering to a great extent; See Section 6.3.1. 6.3 Numerical Examples We start with a number of toy examples to demonstrate the results and then proceed with more realistic and practical examples. 6.3.1 2D System Without Uncertainty Consider the system ẋ = [ 1 1 −1 0 ] x+ [ 1 1 ] u 122 6.3. Numerical Examples with u(t) ∈ U := [−0.15, 0.15] and a constraint set K defined to be a disc of radius 0.5 about the origin. We approximate the viability ker- nel V iab[0,2](K,U) (using |M| = 2 directions and uniform partition with |P | = 80) and synthesize safety-preserving control laws that ensure x(t) ∈ K ∀t ∈ [0, 2]. Since u in qperf can be chosen arbitrarily in U , we simply ap- ply u = 0. The simulation results for a few initial conditions are given in Figures 6.3 and 6.4. Smooth Transition Between Safety and Performance Here we will apply the mildly varying control law described by (6.38)–(6.41) (with α = 0) to obtain less switching frequency for the controller. Fig- ure 6.5 shows the closed-loop trajectories generated by such a policy, while Figure 6.6 shows the corresponding feedback strategies. 6.3.2 2D System With Uncertainty Consider the system ẋ = [ 1 1 −1 0 ] x+ [ 1 1 ] u+ [ 0 1 ] v subject to unknown but bounded disturbance v(t) ∈ V := [−0.1, 0.1], bounded control input u(t) ∈ U := [−0.15, 0.15], and state constraint x(t) ∈ K := {x ∈ R2 | ‖x‖2 ≤ 0.5} for all t ∈ [0, 2]. We approximate the (finite-horizon) discriminating kernel Disc[0,2](K,U ,V) (using |M| = 4 directions and uni- form partition with |P | = 80) and synthesize safety-preserving control laws that ensure x(t) ∈ K ∀t ∈ [0, 2] despite the action of the unknown distur- bance. As before, since u in qperf can be chosen arbitrarily in U , we simply apply u = 0. The simulation results for a few initial conditions are given in Figures 6.7, 6.8, and 6.9, where we have assumed that for each initial condition the disturbance is a random signal with uniform distribution on V, i.e. v(t) ∈ U([−0.1, 0.1]). Safety is preserved despite this (unknown) disturbance. 123 6.3. Numerical Examples −0.6 −0.4 −0.2 0 0.2 0.4 0.6 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 x1 x 2 x0,1 x0,2 x0,3 x0,4 Figure 6.3: Closed-loop trajectories in the phase-plane for four sample ini- tial conditions x0,i for Example 6.3.1. The constraint set K (red/dark) and a piecewise ellipsoidal under-approximation of the viability kernel (green/light) are shown. The level set ap- proximation of the viability kernel is outlined in thick black lines. The feedback law is chosen such that u = 0 when q = qperf and u = usafe when q = qsafe. Safety is preserved over the fi- nite horizon for x0,i ∈ V iab[0,2](K,U). While x0,3 is not in the under-approximation of the viability kernel, the designed control policy still maintains safety (see page 119). Note that the tra- jectories can leave the kernel since it is finite-horizon and hence not invariant; however, they do not violate the constraints. For x0,4 6∈ V iab[0,2](K,U), even though the control is always at its extremum points (see the bottom plot in Figure 6.4), the tra- jectory eventually leaves K. 124 6.3. Numerical Examples 0 0.5 1 1.5 2 −0.1 0 0.1 0 0.5 1 1.5 2 −0.1 0 0.1 0 0.5 1 1.5 2 −0.1 0 0.1 0 0.5 1 1.5 2 −0.1 0 0.1 t u u u u (x0,1) (x0,2) (x0,3) (x0,4) Figure 6.4: The corresponding safety-preserving feedback policy for Exam- ple 6.3.1 for each initial condition x0,i (in Figure 6.3). The dashed lines (red) indicate the hard bounds on the input. 125 6.3. Numerical Examples −0.6 −0.4 −0.2 0 0.2 0.4 0.6 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 x1 x 2 x0,1 x0,2 x0,3 x0,4 Figure 6.5: Closed-loop trajectories with smooth transition (via (6.38)– (6.41) with α = 0) between qperf and qsafe for Example 6.3.1. Safety is maintained. 126 6.3. Numerical Examples 0 0.5 1 1.5 2 −0.1 0 0.1 0 0.5 1 1.5 2 −0.1 0 0.1 0 0.5 1 1.5 2 −0.1 0 0.1 0 0.5 1 1.5 2 −0.1 0 0.1 t u u u u (x0,1) (x0,2) (x0,3) (x0,4) Figure 6.6: The corresponding safety-preserving feedback policy with less switching frequency for Example 6.3.1 for x0,i, i = 1, . . . , 3 (in Figure 6.5). Note that for x0,4 the controller is never in qperf. 127 6.3. Numerical Examples −0.6 −0.4 −0.2 0 0.2 0.4 0.6 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 x1 x 2 x0,1 x0,2 x0,3 Figure 6.7: Closed-loop trajectories in the phase-plane for three sample ini- tial conditions x0,i for Example 6.3.2. The constraint set K (red/dark) and a piecewise ellipsoidal under-approximation of the discriminating kernel (green/light) are shown. The level set approximation of the kernel is outlined in thick black lines. The feedback law is chosen such that u = 0 when q = qperf and u = usafe when q = qsafe. Safety is preserved despite the distur- bance. 128 6.3. Numerical Examples 0 0.5 1 1.5 2 −0.2 0 0.2 0 0.5 1 1.5 2 −0.2 0 0.2 0 0.5 1 1.5 2 −0.2 0 0.2 t u u u (x0,1) (x0,2) (x0,3) Figure 6.8: The corresponding safety-preserving feedback policy for Exam- ple 6.3.2 for each initial condition x0,i. 0 0.5 1 1.5 2 −0.1 0 0.1 0 0.5 1 1.5 2 −0.1 0 0.1 0 0.5 1 1.5 2 −0.1 0 0.1 t v v v (x0,1) (x0,2) (x0,3) Figure 6.9: The simulated disturbance signal for Example 6.3.2 for each ini- tial condition x0,i. 129 6.3. Numerical Examples k10 k12 k13 k21 k31 V1 0.4436 0.1140 0.0419 0.0550 0.0033 16.044 Table 6.1: Model parameters for the patient (11 years old, 35 kg); cf. [1] 6.3.3 3D Control of Anesthesia Consider the problem of safety in control of anesthesia as described in Sec- tion 5.3.2. Instead of the discrete-time Laguerre model of the patient, here we use a three-dimensional continuous-time pharmacokinetic/pharmacody- namic (PKPD) compartmental model whose state variables describe the concentration of propofol in each compartment of the body: ẋ = −(k10 + k12 + k13) k12 k13k21 −k21 0 k31 0 −k31 x+ 1/V10 0 u. The model parameters, taken from [1], are given in Table 6.1. Suppose that the patient is to undergo a 30 min surgery and that the input (propofol infusion rate [mg/min]) is hard bounded above and below such that u(t) ∈ U := [0, 20]. To keep the plasma drug concentration within therapeutic range, we require x(t) ∈ K := E 3.55 5 , 6.25 0 00 25 0 0 0 25 ∀t ∈ [0, 30]. (6.43) We approximate the viability kernel V iab[0,30](K,U) via Algorithm 5.1 using |M| = 15 directions and uniform partition with |P | = 600. Consider the initial condition x0 = [ 4 5 8 ]T and assume that no drug is being ad- ministered, i.e. u = min(U) = 0. When the safety-preserving controller is not engaged, the therapeutic constraint K is eventually violated (Figure 6.10) putting the patient at risk of intra- and post-operative complications. When the safety-preserving controller is engaged (using the same input u = 0 in qperf and enforcing the optimal safety control law in qsafe) the constraint K is satisfied for the entire duration of surgery (Figure 6.11). To ensure 130 6.3. Numerical Examples 0 2 4 6 0 5 10 0 2 4 6 0 5 10 0 5 10 0 5 10 x1 x 2 x1 x 3 x 3 x2 Figure 6.10: 2D projections of the closed-loop trajectories with u = 0 and no safety control for Example 6.3.3. The initial condition x0 is marked by ‘×’. The constraint set K (red/dark) and a piecewise ellipsoidal under-approximation of the viability ker- nel (green/light) are also shown. The constraint set K (safety of the patient) is violated. a slowly varying safety-preserving infusion policy the modified strategy in (6.38)–(6.41) (with α = 0) is employed: Figure 6.12 shows this feedback policy, while Figures 6.13–6.16 discuss various behaviors and characteristics of the hybrid controller using such a strategy. In comparison, a direct appli- cation of the safety-preserving controller without the modified strategy in (6.38)–(6.41), while still capable of maintaining safety, results in significant chattering and an aggressive infusion policy as shown in Figures 6.17–6.21. 6.3.4 6D Flight Envelope Protection Consider the problem of aerodynamic flight envelope protection for NASA’s Highly Maneuverable Aircraft Technology (HiMAT) [105]. The longitudinal dynamics of the HiMAT aircraft trimmed at 25 kft and 0.9 Mach are given 131 6.3. Numerical Examples 0 2 4 6 0 5 10 0 2 4 6 0 5 10 0 5 10 0 5 10 x2 x 3 x 2 x1 x1 x 3 Figure 6.11: 2D projections of the closed-loop trajectories with u = usafe in qsafe and uperf = 0 in qperf using the modified policy (6.38)– (6.41) with α = 0 for Example 6.3.3. The initial condition x0 is marked by ‘×’. The constraint set K (red/dark) and a piecewise ellipsoidal under-approximation of the viability ker- nel (green/light) are also shown. Safety is preserved. 0 5 10 15 20 25 30 0 5 10 15 20 u t Figure 6.12: The safety-preserving feedback policy via (6.38)–(6.41) with α = 0 for Figure 6.11 for Example 6.3.3. The bounds on the input are shown as dashed lines. 132 6.3. Numerical Examples 0 5 10 15 20 25 30 0 20 40 60 80 100 % u pe rf in u t Figure 6.13: uperf(t) as a percentage of u(t) in qperf for the modified policy (6.38)–(6.41) with α = 0 for Example 6.3.3. The safety com- ponent of the input becomes more dominant as the trajectory gets closer to the boundary of the ellipsoid that is the domain of qperf as shown in Figure 6.16. A value of 0% would corre- spond to when the automaton has switched to qsafe, while a value of 100% would correspond to when x(t) is at the center of the domain of qperf. Increasing α would further emphasize the uperf component, at the cost of a more aggressive control policy and possibly more frequent switchings between the two modes of the automaton. 0 5 10 15 20 25 30 perf safe a ct iv e m od e t Figure 6.14: The active mode of the hybrid automaton H at time t using the modified policy (6.38)–(6.41) with α = 0 for Example 6.3.3. The safety component usafe of the input in qperf seems to be sufficient to maintain safety in this case, and thus the controller does not switch to qsafe. 133 6.3. Numerical Examples 0 5 10 15 20 25 30 5 10 15 γ t Figure 6.15: The index of the ellipsoid being used by the automaton (i.e. the internal input γ ∈ Σi of H) at time t when using the modified policy (6.38)–(6.41) with α = 0 for Example 6.3.3. The ellipsoid does not change. 0 5 10 15 20 25 30 0 0.2 0.4 0.6 0.8 1 φ t Figure 6.16: The location of x(t) within the domain of qperf, i.e. φ [γ](x(t), t), when using the modified policy (6.38)–(6.41) with α = 0 for Example 6.3.3. The boundary of the corresponding ellipsoid is marked by the dashed line at 1. The center of the ellipsoid is at 0. 134 6.3. Numerical Examples 0 2 4 6 0 5 10 0 2 4 6 0 5 10 0 5 10 0 5 10 x2 x 2 x1x1 x 3 x 3 Figure 6.17: 2D projections of the closed-loop trajectories with u = uperf = 0 in qperf and u = usafe in qsafe for Example 6.3.3. While safety is still preserved, the corresponding feedback policy chatters significantly (Figure 6.18) as compared to the slowly varying policy in Figure 6.12. 0 5 10 15 20 25 30 0 5 10 15 20 u t Figure 6.18: The safety-preserving feedback policy for Figure 6.17 for Ex- ample 6.3.3. The bounds on the input are shown as dashed lines. The controller chatters between qperf and qsafe as shown in Figure 6.19. 135 6.3. Numerical Examples 0 5 10 15 20 25 30 perf safe a ct iv e m od e t Figure 6.19: The active mode of the hybrid automaton H at time t for Ex- ample 6.3.3 when the modified policy (6.38)–(6.41) is not em- ployed. 0 5 10 15 20 25 30 5 10 15 γ t Figure 6.20: The index of the ellipsoid being used by the automaton (i.e. the internal input γ ∈ Σi of H) at time t for Example 6.3.3. Note that the ability to switch ellipsoids in H is not the cause of the chatter (as evidenced by the fact that the ellipsoid appears to only switch twice). 136 6.3. Numerical Examples 0 5 10 15 20 25 30 0 0.2 0.4 0.6 0.8 1 φ t Figure 6.21: The location of x(t) within the domains of qperf and qsafe, i.e. φ[γ](x(t), t), for Example 6.3.3. For φ[γ](x(t), t) < 1 the active mode is qperf (whose domain is an ellipsoid whose index is shown in Figure 6.20). On the other hand, for φ[γ](x(t), t) ≥ 1 the active mode is qsafe. At time t = 1.71 the automaton takes a transition on (qperf, qperf) ∈ E while resetting γ to correspond to a new direction vector (ellipsoid). This ensures that the active mode remains qperf. The change of ellipsoid causes a change in the definition of φ and hence a discontinuous jump in the plot. At t = 4.1 the automaton is forced to switch to qsafe to maintain safety, but progresses by chattering between qsafe and qperf as indicated by the values of φ [γ](x(t), t). 137 6.3. Numerical Examples by ẋ = Ax+Bu with A = −0.0226 −36.6170 −18.8970 −32.0900 3.2509 −0.7626 0.0001 −1.8997 0.9831 −0.0007 −0.1708 −0.0050 0.0123 11.7200 −2.6316 0.0009 −31.6040 22.3960 0 0 1.0000 0 0 0 0 0 0 0 −30.0000 0 0 0 0 0 0 −30.0000 , B = [ 0 0 0 0 30 0 0 0 0 0 0 30 ]T and state vector x = [ α̇ α θ̇ θ xδe xδc ]T in which the first four states represent angle of attack and attitude angle and their rates of change, and the last two states represent elevon and canard control actuator dynamics. The control input u = [ δe δc ]T is comprised of elevon and canard actuators. For all t ∈ [0, 1] we assume that u(t) ∈ U , a disc of radius 0.5236 rad (≈ 30◦) about the origin in R2, and require x(t) ∈ K := {x ∈ R6 | ‖x‖2 ≤ 5}.1 Suppose that a pre-designed Linear Quadratic Regulator (LQR) con- troller is to be used that satisfies the performance functional J(u) = ∫ ∞ 0 (xTQx+ uTRu)dt (6.44) with Q = I6 and R = 10 2 × I2 subject to the system dynamics. The corresponding state-feedback control law that minimizes J(u) is ulqr := −Kx with K = [ 0.0823 −0.8209 −0.3457 −0.5247 0.3119 −0.2107 −0.0552 0.5619 0.2367 0.3588 −0.2107 0.1497 ] . Since the control authority is constrained by the ellipsoid U =: E(ω,Ω), 1These assumptions are purely for the algorithmic convenience of dealing with ellip- soids. The input constraint set, for instance, may be better described by a rectangle in R2 as the actuators are in fact dynamically independent, and the state constraint set is entirely fictitious. 138 6.4. Summary and Further Discussions applying the above LQR control law may result in saturation such that u = sat(ulqr) is effectively applied to the system. Here the saturation function sat : R2 → U is determined by the support vector of the set U in the direction of ulqr/‖ulqr‖ and is defined as sat(ulqr) := ulqr if ulqr ∈ U ;ω + Ωulqr〈ulqr, Ωulqr〉−1/2 if ulqr 6∈ U . (6.45) We approximate the viability kernel V iab[0,1](K,U) via Algorithm 5.1 using |M| = 15 directions and uniform partition with |P | = 300. Consider the initial condition x0 = [ −1.7064 1.7769 −1.8770 −1.1272 1.5994 1.7680 ]T . Applying the LQR controller without a safety-preserving strategy in place results in violation of the aerodynamic envelope K as shown in Figure 6.22 due to unaccounted actuator saturations. On the other hand, when the safety-preserving controller is employed (with u = sat(ulqr) in mode qperf and the optimal safety control law in qsafe) the flight envelope is protected over the horizon [0, 1] despite the unaccounted actuator saturations in the LQR action (Figure 6.23). 6.4 Summary and Further Discussions In this chapter we first extended the results of Chapter 5 to the differential games setting in which the system is subject to unknown but bounded dis- turbance/uncertainty. Consequently, the discriminating kernel is expressed in terms of robust maximal reachable sets. Owing to this new connec- tion, scalable and efficient Lagrangian methods can now be used to cor- rectly approximate the discriminating kernel. We also presented an exten- sion of the piecewise ellipsoidal algorithm from Section 5.2.1 that facilitates the under-approximation of the discriminating kernel for high-dimensional LTI systems. Based on this algorithm (and its underlying ellipsoidal tech- 139 6.4. Summary and Further Discussions -20 -10 0 10 -5 0 5 -20 -10 0 10 -5 0 5 -20 -10 0 10 -5 0 5 -5 0 5 -5 0 5 -5 0 5 -5 0 5 -5 0 5 -5 0 5 x1 x1 x1 x 2 x2 x2 x 3 x 3 x3 x 4 x 4 x 4 Figure 6.22: 2D projections of the closed-loop trajectories with saturated LQR and no safety control for the first four states for Ex- ample 6.3.4. The initial condition x0 is marked by ‘×’. The constraint set K (red/dark) and a piecewise ellipsoidal under- approximation of the viability kernel (green/light) are also shown. The constraint set K is violated due to unaccounted actuator saturations. 140 6.4. Summary and Further Discussions -5 0 5 -5 0 5 -5 0 5 -5 0 5 -5 0 5 -5 0 5 -5 0 5 -5 0 5 -5 0 5 -5 0 5 -5 0 5 -5 0 5 x 2 x 4 x 4 x 3 x 3 x 4 x2 x3 x1 x2 x1 x1 Figure 6.23: 2D projections of the closed-loop trajectories with safety- preserving controller for the first four states for Example 6.3.4. The initial condition x0 is marked by ‘×’. The feedback law is chosen such that u = sat(ulqr) when q = qperf and u = usafe when q = qsafe. The constraint set K (red/dark) and a piece- wise ellipsoidal under-approximation of the viability kernel (green/light) are shown. Safety is preserved in spite of the actuator saturations in the LQR controller. 141 6.4. Summary and Further Discussions niques for reachability [72]), we then proposed a scalable, non-restrictive, safety-preserving, hybrid state-feedback control strategy for continuous-time LTI systems that ensures that the state constraint is not violated despite bounded control authority and, if present, unknown disturbances. We demon- strated the results on several examples including a realistic problem of safety in anesthesia and a 6D problem of aerodynamic flight envelope protection. 142 Chapter 7 Conclusions and Future Work We considered the problem of guaranteed safety and constraint satisfaction in high-dimensional, safety-critical, controlled dynamic systems. Reacha- bility analysis and viability theory provide an appropriate framework for set-valued analysis of constrained dynamical systems. To guarantee safety of such systems and to synthesize controllers that are capable of preserving this safety despite bounded control authority (and possibly disturbances or uncertainties), the computation of the minimal reachable tube or by duality, the viability kernel is required. Historically, the algorithms that approxi- mate these sets—known as Eulerian methods—are based on gridding the state space. While powerful and versatile, their computational complexity increases exponentially with the dimension of the state. This renders them impractical for systems of dimensions higher than three or four. We presented two separate approaches for reduction of complexity in computing the minimal reachable tube or the viability kernel for higher- dimensional systems. The first approach, based on structure decomposi- tion, aims to facilitate the use of Eulerian methods on higher-dimensional, continuous-time, continuously-valued LTI systems (and by extension, hy- brid systems with LTI continuous dynamics). It does so by constructing an appropriate similarity transformation that not only results in decoupling (or weak unidirectional coupling) of the dynamics, but also yields disjoint input. This imposed structure is then exploited for decomposition of the system for the purpose of computing the minimal reachable tube and the viability kernel in lower dimensions. A number of algorithms are then presented that 143 7.1. Summary of Contributions enable a sound approximation of these constructs in lower-dimensional sub- spaces. It is shown that the reverse transformation of the intersection of the back-projection of these resulting sets correctly approximates the actual con- structs. Within the framework of structure decomposition, we proposed two techniques: the Schur-based decomposition and the Riccati-based decompo- sition, each with its own merits. While the Schur-based decomposition is quite generic and applicable to most systems, the Riccati-based decomposi- tion may yield less conservative reachability computations for two-time-scale or ill-conditioned systems as was shown with an example in Chapter 4. The second complexity reduction approach, based on set-theoretic meth- ods, draws a connection between the viability kernel and the maximal reach- able sets for continuous- and discrete-time systems. Since the maximal reachable sets can be computed using efficient and scalable techniques— known as Lagrangian methods—that employ compact set representations and follow the flow of the dynamics, the viability kernel can now be under- approximated with polynomial complexity. Using the well-established el- lipsoidal techniques for maximal reachability we then proposed a scalable algorithm that facilitates the computation of the viability kernel for high- dimensional LTI systems. This approach also enabled us to formulate a scalable safety-preserving static feedback control strategy. We also provided extensions of this approach to systems with unknown but bounded distur- bances or uncertainties. We demonstrated our techniques on several examples (up to 8D) includ- ing a problem of safety in control of anesthesia and flight envelope protection for longitudinal aircraft dynamics. We mention, however, that the second approach is much more scalable: The algorithm can likely be applied to systems with several dozens of state variables. 7.1 Summary of Contributions We summarize our contributions as follows: • We proposed two structure decomposition techniques that enable the 144 7.2. Future Research Directions computation of the minimal reachable tube and the viability kernel for higher-dimensional LTI systems using Eulerian methods. The decom- position techniques that are available in the literature, while generally applicable to nonlinear systems, assume a certain structure on the system that can be exploited. Our techniques, on the other hand, are designed to impose this structure. • We presented a novel connection between the viability kernel and max- imal reachable sets. Owing to this connection efficient and scalable Lagrangian methods can now be used to approximate the viability kernel for high-dimensional systems. Our piecewise ellipsoidal algo- rithm which was proposed based on this new connection using the ellipsoidal techniques for maximal reachability is capable of comput- ing a guaranteed under-approximation of the viability (discriminating) kernel for high-dimensional (uncertain) LTI systems. • We showed that the presented connection can also be employed to synthesize safety-preserving control laws. We then proposed a scal- able safety-preserving control strategy (again based on the ellipsoidal techniques for maximal reachability and the corresponding optimal control laws) that ensures safety of high-dimensional safety-critical, possibly uncertain LTI systems. 7.2 Future Research Directions There are several possibilities for future work and further developments. Structure Decomposition For both of our decomposition techniques—i.e. Schur-based decomposition (Chapters 3) and Riccati-based decomposition (Chapter 4)—future work includes efforts to reduce potential conservatism in the over-approximation of the minimal reachable tube (or the under-approximation of the viability kernel). One direction is to incorporate the geometric information about the 145 7.2. Future Research Directions shape of the target set into the decomposition process so that the projection of the set onto the subspaces of the transformed coordinates does not result in excessive loss of detail. A second direction is an alternative transformation that produces subsystems that may both be manipulated to some degree while still preserving the disjoint property of the input. Set-Theoretic Methods While the presented piecewise ellipsoidal algorithm in Chapter 5 has proven to be effective and efficient, it may be subject to excessive conservatism par- ticularly for large time horizons. Quantifying the accuracy of the algorithm is an important future research direction. We are also currently developing alternative approaches that yield a more accurate under-approximation of the viability kernel while still preserving the scalability property. In fact our hope is that the presented connection between the viability kernel and maximal reachable sets encourages the development of scalable and accurate algorithms for the computation of the viability kernel for nonlinear systems. Finally, as our algorithm is already heavily based on intersections, a natu- ral (and straightforward) extension would be to consider hybrid dynamical systems. Safety-Preserving Control Synthesis Any Lagrangian technique that can accommodate the synthesis of maximal reachability control laws (and can handle adversarial inputs for the discrim- inating kernel case) may be used to formulate a safety-preserving controller based on the framework described in Chapter 6. Our proposed control al- gorithm in that chapter is at early stages of development. Many future research directions are possible: • We have assumed that safety is only to be preserved over the given finite interval [0, τ ]—beyond this point the trajectories may leave the constraint set regardless of what the control input does. We are cur- rently developing a variation of the hybrid controller (6.18) that em- 146 7.2. Future Research Directions ploys a pseudo-time variable with varying rates of change and can potentially ensure safety over a horizon larger than τ . • Suppose that the infinite-horizon discriminating kernel is nonempty, and Disc[0,τ−tk](K,U ,V) ≡ DiscR+(K,U ,V) for some k ∈ {0, . . . , |P |}. Due to propagation of the approximation error, the set generated by the piecewise ellipsoidal algorithm is not in general guaranteed to con- verge to the infinite-horizon kernel. However, if the algorithm does converge for at least one terminal direction, then the current formula- tion of the hybrid controller can be used to synthesize infinite-horizon safety-preserving control laws. We intend to prove this point in the future. • Another avenue is the extension of the safety-preserving controller to the discrete-time case. While this is certainly possible, the discussions become more involved in the presence of disturbance/uncertainty as the Isaac’s condition no longer holds. • The presented safety-preserving controller is non-restrictive/permissive. That is, the optimal safety control law must be applied only when safety is at stake. Otherwise, any desired/performance-satisfying con- trol law in U can be chosen. A future research direction would be to determine the “distance” of this controller (in an ordered space of all safety-preserving controllers) to the least restrictive controller as defined by [81]. 147 Bibliography [1] A. Absalom and G. Kenny. Paedfusor pharmacokinetic data set. British Journal of Anaesthesia, 95(1):110, 2005. [2] E. Asarin, O. Bournez, T. Dang, O. Maler, and A. Pnueli. Effective synthesis of switching controllers for linear systems. Proceedings of the IEEE, 88(7):1011–1025, 2000. [3] E. Asarin and T. Dang. Abstraction by projection and application to multi-affine systems. In Hybrid Systems: Computation and Control, LNCS 2993, pages 129–132. Springer-Verlag, 2004. [4] E. Asarin, T. Dang, G. Frehse, A. Girard, C. Le Guernic, and O. Maler. Recent progress in continuous and hybrid reachability analy- sis. In Proc. IEEE International Symposium on Computer-Aided Con- trol Systems Design, Munich, Germany, October 2006. [5] J.-P. Aubin. Viability Theory. Systems and Control: Foundations and Applications. Birkhäuser, Boston, MA, 1991. [6] J.-P. Aubin. Viability kernels and capture basins of sets under differ- ential inclusions. In Proc. IEEE Conference on Decision and Control, pages 4605–4610, Las Vegas, NV, Dec 2002. [7] J.-P. Aubin, A. M. Bayen, and P. Saint-Pierre. Viability Theory: New Directions. Springer Verlag, 2nd edition, 2011. [8] A. M. Bayen, I. M. Mitchell, M. Oishi, and C. J. Tomlin. Aircraft au- tolander safety analysis through optimal control-based reach set com- putation. Journal of Guidance, Control, and Dynamics, 30(1):68–77, 2007. [9] C. Béné, L. Doyen, and D. Gabay. A viability analysis for a bio- economic model. Ecological Economics, 36(3):385–396, 2001. 148 Bibliography [10] S. Bibian, C. Ries, M. Huzmeman, and G. A. Dumont. Introduction to automated drug delivery in clinical anesthesia. European Journal of Control, 11(6):535–557, Dec 2005. [11] F. Blanchini. Set invariance in control. Automatica, 35(11):1747–1767, 1999. [12] F. Blanchini and S. Miani. Set-Theoretic Methods in Control. Springer, 2008. [13] F. Borrelli, C. Del Vecchio, and A. Parisio. Robust invariant sets for constrained storage systems. Automatica, 45(12):2930–2936, 2009. [14] G. Bouligand. Introducion a la geometrie inxnitesimale directe. Paris, Bourgin: Gauthiers-Villars, 1932. [15] S. P. Boyd and L. Vandenberghe. Convex optimization. Cambridge University Press, 2004. [16] A. E. Bryson. Control of Spacecraft and Aircraft. Princeton Univ. Press, 1994. [17] J. Camara, A. Girard, and G. Gossler. Safety controller synthesis for switched systems using multi-scale symbolic models. In Proc. Con- ference on Decision and Control and European Control Conference, pages 520–525, Orlando, FL, 2011. [18] P. Cardaliaguet, M. Quincampoix, and P. Saint-Pierre. Set-valued numerical analysis for optimal control and differential games. In M. Bardi, T. Raghavan, and T. Parthasarathy, editors, Stochastic and Differential Games: Theory and Numerical Methods, number 4 in An- nals of the International Society of Dynamic Games, pages 177–247, Boston, MA, 1999. Birkhäuser. [19] K. W. Chang. Singular perturbations of a general boundary value problem. SIAM Journal on Mathematical Analysis, 3:520–526, 1972. [20] A. Chutinan and B. H. Krogh. Computing polyhedral approximations to flow pipes for dynamic systems. In Proc. IEEE Conference on Decision and Control, pages 2089–2094, Tampa, FL, December 1998. [21] P.-A. Coquelin, S. Martin, and R. Munos. A dynamic programming approach to viability problems. In Proc. IEEE Symposium on Approx- imate Dynamic Programming and Reinforcement Learning (ADPRL 2007), pages 178–184, 2007. 149 Bibliography [22] T. Dang and O. Maler. Reachability analysis via face lifting. In Hy- brid Systems: Computation and Control, LNCS 1386, pages 96–109. Springer, 1998. [23] A. N. Daryin, A. B. Kurzhanski, and I. V. Vostrikov. Reachability ap- proaches and ellipsoidal techniques for closed-loop control of oscillating systems under uncertainty. In Proc. IEEE Conference on Decision and Control, pages 6385–6390, San Diego, CA, 2006. [24] G. Deffuant, L. Chapel, and S. Martin. Approximating viability ker- nels with support vector machines. IEEE Transactions on Automatic Control, 52(5):933–937, 2007. [25] D. Del Vecchio, M. Malisoff, and R. Verma. A separation principle for a class of hybrid automata on a partial order. In Proc. American Control Conference, pages 3638–3643, 2009. [26] G. Dumont, A. Martinez, and J. Ansermino. Robust control of depth of anesthesia. International Journal of Adaptive Control and Signal Processing, 23:435–454, 2009. [27] L. Evans and P. Souganidis. Differential games and representation formulas for solutions of Hamilton-Jaconbi-Isaacs equations. Indiana University Mathematics Journal, 33(5):773–797, 1984. [28] F. Fadaie and M. E. Broucke. On the least restrictive control for collision avoidance of two unicycles. International Journal of Robust and Nonlinear Control, 16(12):553–574, 2006. [29] A. F. Filippov and F. M. Arscott. Differential Equations With Discon- tinuous Righthand Sides. Mathematics and Its Applications. Kluwer Academic Publishers, 1988. [30] 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 G. Gopalakrishnan and S. Qadeer, editors, Proc. 23rd International Conference on Computer Aided Ver- ification (CAV), pages 1–16. Springer, 2011. [31] G. Freiling. A survey of nonsymmetric Riccati equations. Linear Algebra and its Applications, 351:243–270, 2002. 150 Bibliography [32] Z. Gajic and I. Borno. General transformation for block diagonal- ization of weakly coupled linear systems composed of N-subsystems. IEEE Transactions on Circuits and Systems—Part I: Fundamental Theory and Applications, 47(6):909–912, 2000. [33] Y. Gao, J. Lygeros, and M. Quincampoix. The reachability problem for uncertain hybrid systems revisited: a viability theory perspective. In J. Hespanha and A. Tiwari, editors, Hybrid Systems: Computation and Control, LNCS 3927, pages 242–256, Berlin Heidelberg, 2006. Springer-Verlag. [34] T. Gilhuly. Modeling and control of neuromuscular blockade. PhD thesis, University of British Columbia, Vancouver, Canada, 2007. [35] T. Gilhuly, G. Dumont, and B. Macleod. Modelling for computer controlled neuromuscular blockade. In IEEE Engineering in Medicine and Biology Magazine, pages 26–29, 2005. [36] A. Girard. Reachability of uncertain linear systems using zonotopes. In M. Morari, L. Thiele, and F. Rossi, editors, Hybrid Systems: Com- putation and Control, LNCS 3414, pages 291–305. Springer, 2005. [37] A. Girard. Synthesis using approximately bisimilar abstractions: state-feedback controllers for safety specifications. In Hybrid Systems: Computation and Control, Stockholm, Sweden, 2010. [38] A. Girard. Controller synthesis for safety and reachability via approx- imate bisimulation. Automatica, 48(5):947–953, 2012. [39] A. Girard and C. Le Guernic. Efficient reachability analysis for linear systems using support functions. In IFAC World Congress, Seoul, Korea, July 2008. [40] A. Girard, C. Le Guernic, and O. Maler. Efficient computation of reachable sets of linear time-invariant systems with inputs. In J. Hes- panha and A. Tiwari, editors, Hybrid Systems: Computation and Con- trol, LNCS 3927, pages 257–271. Springer-Verlag, 2006. [41] A. Girard and G. J. Pappas. Approximate bisimulation relations for constrained linear systems. Automatica, 43(8):1307–1317, 2007. [42] A. Girard and G. J. Pappas. Approximation metrics for discrete and continuous systems. IEEE Transactions on Automatic Control, 52(5):782–798, 2007. 151 Bibliography [43] G. H. Golub and C. F. Van Loan. Matrix Computations. Johns Hop- kins Univ. Press, 1996. [44] M. R. Greenstreet. Verifying safety properties of differential equations. In Proc. Conference on Computer Aided Verification, pages 277–287, New Brunswick, NJ, July 1996. [45] J. Groß. Explicit solutions to the matrix inverse problem AX = B. Linear Algebra and its Applications, 289:131–134, 1999. [46] M. R. Hafner, D. Cunningham, L. Caminiti, and D. Del Vecchio. Auto- mated vehicle-to-vehicle collision avoidance at intersections. In Proc. of World Congress on Intelligent Transport Systems, Orlando, FL, October 2011. [47] Z. Han and B. H. Krogh. Reachability analysis of hybrid control sys- tems using reduced-order models. In Proc. American Control Confer- ence, pages 1183–1189, Boston, MA, June 2004. [48] Z. Han and B. H. Krogh. Reachability analysis for affine systems using -decomposition. In Proc. IEEE Conference on Decision and Control, and European Control Conference, pages 6984–6990, Seville, Spain, December 2005. [49] Z. Han and B. H. Krogh. Reachability analysis of large-scale affine systems using low-dimensional polytopes. In J Hespanha and A. Ti- wari, editors, Hybrid Systems: Computation and Control, LNCS 3927, pages 287–301, Berlin, Germany, 2006. Springer-Verlag. [50] Z. Han and B. H. Krogh. Reachability analysis of nonlinear systems using trajectory piecewise linearized models. In Proc. American Con- trol Conference, pages 1505–1510, Minneapolis, MN, 2006. [51] P. A. Ioannou and J. Sun. Robust Adaptive Control. Prentice Hall, Englewood Cliffs, NJ, 1996. [52] C. Ionescu, R. De Keyser, B. Torrico, T. De Smet, M. Struys, and J. Normey-Rico. Robust predictive control strategy applied for propo- fol dosing using BIS as a controlled variable during anesthesia. IEEE Transactions on Biomedical Engineering, 55(9):2161–2170, 2008. [53] K. J. Johansson, M. Egerstedt, J. Lygeros, and S. Sastry. On the regularization of Zeno hybrid automata. Systems and Control Letters, 38:141–150, 1999. 152 Bibliography [54] A. Kanade, R. Alur, F. Ivančić, and S. Ramesh. Generating and analyzing symbolic traces of Simulink/Stateflow models. In A. Bouaj- jani and O. Maler, editors, Computer Aided Verification, LNCS 5643, pages 430–445. Springer Berlin Heidelberg, 2009. [55] S. Kaynama, J. Maidens, M. Oishi, I. M. Mitchell, and G. A. Dumont. Computing the viability kernel using maximal reachable sets. In Hy- brid Systems: Computation and Control, pages 55–63, Beijing, China, 2012. [56] S. Kaynama and M. Oishi. A modified Riccati transformation for com- plexity reduction in reachability analysis of linear time-invariant sys- tems. IEEE Transactions on Automatic Control. (accepted; preprint available at www.ece.ubc.ca/∼kaynama). [57] S. Kaynama and M. Oishi. Schur-based decomposition for reachability analysis of linear time-invariant systems. In Proc. IEEE Conference on Decision and Control, pages 69–74, Shanghai, China, December 2009. [58] S. Kaynama and M. Oishi. Complexity reduction through a Schur- based decomposition for reachability analysis of linear time-invariant systems. International Journal of Control, 84(1):165–179, 2011. [59] S. Kaynama, M. Oishi, I. M. Mitchell, and G. A. Dumont. The con- tinual reachability set and its computation using maximal reachability techniques. In Proc. IEEE Conference on Decision and Control, and European Control Conference, pages 6110–6115, Orlando, FL, 2011. [60] S. Kaynama, M. Oishi, I. M. Mitchell, and G. A. Dumont. Fixed- complexity piecewise ellipsoidal representation of the continual reach- ability set based on ellipsoidal techniques. In Proc. American Control Conference, Montreal, QC, 2012. (to appear). [61] E. C. Kerrigan. Robust Constraint Satisfaction: Invariant Sets and Predictive Control. PhD thesis, University of Cambridge, 2000. [62] P. V. Kokotović. A Riccati equation for block-diagonalization of ill-conditioned systems. IEEE Transactions on Automatic Control, 20(6):812–814, 1975. [63] P. V. Kokotović, H. K. Khalil, and J. O’Reilly. Singular Perturbation Methods in Control: Analysis and Design. SIAM, 1999. 153 Bibliography [64] E. K. Kostousova. Control synthesis via parallelotopes: optimiza- tion and parallel computations. Optimization methods and software, 14(4):267–310, 2001. [65] B. H. Krogh and O. Stursberg. Efficient representation and computa- tion of reachable sets for hybrid systems. In O. Maler and A. Pnueli, editors, Hybrid Systems: Computation and Control, LNCS 2623, pages 482–497, Berlin, Germany, 2003. Springer-Verlag. [66] W. Kühn. Rigorously computed orbits of dynamical systems without the wrapping effect. Computing, 61(1):47–67, 1998. [67] A. B. Kurzhanski and T.F. Filippova. On the description of the set of viable trajectories of a differential inclusion. Sov. Math. Doklady, 34, 1987. [68] A. B. Kurzhanski, I. M. Mitchell, and P. Varaiya. Optimization tech- niques for state-constrained control and obstacle problems. Journal of Optimization Theory and Applications, 128(3):499–521, 2006. [69] A. B. Kurzhanski and I. Vályi. Ellipsoidal Calculus for Estimation and Control. Birkhäuser, Boston, MA, 1996. [70] A. B. Kurzhanski and P. Varaiya. Ellipsoidal techniques for reacha- bility analysis. In N. Lynch and B. Krogh, editors, Hybrid Systems: Computation and Control, LNCS 1790, pages 202–214, Berlin Heidel- berg, 2000. Springer-Verlag. [71] A. B. Kurzhanski and P. Varaiya. Ellipsoidal techniques for reacha- bility analysis: internal approximation. Systems & Control Letters, 41:201–211, 2000. [72] A. B. Kurzhanski and P. Varaiya. On reachability under uncertainty. SIAM Journal on Control and Optimization, 41(1):181–216, 2002. [73] A. A. Kurzhanskiy and P. Varaiya. Ellipsoidal Toolbox (ET). In Proc. IEEE Conference on Decision and Control, pages 1498–1503, San Diego, CA, December 2006. [74] A. A. Kurzhanskiy and P. Varaiya. Ellipsoidal techniques for reacha- bility analysis of discrete-time linear systems. IEEE Transactions on Automatic Control, 52(1):26–38, 2007. 154 Bibliography [75] A. A. Kurzhanskiy and P. Varaiya. Computation of reach sets for dynamical systems. In Chapter for Control Handbook. 2 edition, 2008. [76] M. Kvasnica, P. Grieder, M. Baotić, and M. Morari. Multi-Parametric Toolbox (MPT). In R. Alur and G. J. Pappas, editors, Hybrid Sys- tems: Computation and Control, LNCS 2993, pages 448–462, Berlin, Germany, 2004. Springer. [77] C. Le Guernic. Reachability analysis of hybrid systems with linear continuous dynamics. PhD thesis, Université Grenoble 1 – Joseph Fourier, 2009. [78] C. Le Guernic and A. Girard. Reachability analysis of linear sys- tems using support functions. Nonlinear Analysis: Hybrid Systems, 4(2):250–262, 2010. [79] J. Lygeros. On reachability and minimum cost optimal control. Au- tomatica, 40(6):917–927, June 2004. [80] J. Lygeros, D. N. Godbole, and S. Sastry. Verified hybrid controllers for automated vehicles. IEEE Transactions on Automatic Control, 43(4):522–539, Apr 1998. [81] J. Lygeros, C. J. Tomlin, and S. Sastry. Controllers for reachability specifications for hybrid systems. Automatica, 35:349–370, 1999. [82] M. S. Mahmoud and M. G. Singh. Large Scale Systems Modelling. Pergamon Press, 1981. [83] K. Margellos and J. Lygeros. Air traffic management with target windows: An approach using reachability. In Proc. IEEE Conference on Decision and Control, pages 145–150, Shanghai, China, Dec 2009. [84] D. Q. Mayne, J. B. Rawlings, C. V. Rao, and P. O. Scokaert. Con- strained model predictive control: Stability and optimality. Automat- ica, 36:789–814, 2000. [85] T. Mendonca, J. Lemos, H. Magalhaes, P. Rocha, and S. Esteves. Drug delivery for neuromuscular blockade with supervised multimodel adaptive control. IEEE Transactions on Control Systems Technology, 17(6):1237–1244, November 2009. [86] I. M. Mitchell. Comparing forward and backward reachability as tools for safety analysis. In A. Bemporad, A. Bicchi, and G. Buttazzo, 155 Bibliography editors, Hybrid Systems: Computation and Control, LNCS 4416, pages 428–443, Berlin Heidelberg, 2007. Springer-Verlag. [87] I. M. Mitchell. A toolbox of level set methods. Technical report, UBC Department of Computer Science, TR-2007-11, June 2007. [88] I. M. Mitchell. Scalable calculation of reach sets and tubes for nonlin- ear systems with terminal integrators: a mixed implicit explicit for- mulation. In Proc. Hybrid Systems: Computation and Control, pages 103–112, Chicago, IL, 2011. ACM. [89] 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. [90] 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. [91] I. M. Mitchell and C. J. Tomlin. Overapproximating reachable sets by Hamilton-Jacobi projections. Journal of Scientific Computing, 19(1– 3):323–346, 2003. [92] A. S. Morse. Supervisory control of families of linear set-point controllers—part 1: exact matching. IEEE Transactions on Auto- matic Control, 41:1413–1431, 1996. [93] T. S. Motzkin, H. Raiffa, G. L. Thompson, and R. M. Thrall. The double description method. In H. W. Kuhn and A. W. Tucker, editors, Contributions to the Theory of Games, volume II, pages 51–73, 1953. [94] G. Nenninger, G. Frehse, and V. Krebs. Reachability analysis and control of a special class of hybrid systems. In Modelling, Analysis, and Design of Hybrid Systems, LNCIS 279, pages 173–192. Springer Berlin Heidelberg, 2002. [95] J. Nocedal and S. J. Wright. Numerical Optimization. Springer, New York, NY, 1999. 156 Bibliography [96] M. Oishi, C. J. Tomlin, V. Gopal, and D. N. Godbole. Addressing multiobjective control: Safety and performance through constrained optimization. In Hybrid Systems: Computation and Control, LNCS 2034, pages 459–472. Springer-Verlag, 2001. [97] P. Oliveira, J. P. Hespanha, J. M. Lemos, and T. Mendonça. Super- vised multi-model adaptive control of neuromuscular blockade with off-set compensation. In Proc. European Control Conference, 2009. [98] D. Panagou, K. Margellos, S. Summers, J. Lygeros, and K. J. Kyri- akopoulos. A viability approach for the stabilization of an underac- tuated underwater vehicle in the presence of current disturbances. In Proc. IEEE Conference on Decision and Control, pages 8612–8617, December 2009. [99] S. Prajna and A. Jadbabaie. Safety verification of hybrid systems using barrier certificates. In R. Alur and G. Pappas, editors, Hybrid Systems: Computation and Control, volume LNCS 2993, pages 477– 492, 2004. [100] S. V. Raković. Set theoretic methods in model predictive control. Nonlinear Model Predictive Control, pages 41–54, 2009. [101] S. V. Raković and D. Q. Mayne. Set robust control invariance for linear discrete time systems. pages 975–980, 2005. [102] A. Rantzer and S. Prajna. On analysis and synthesis of safe control laws. In Proc. the Allerton Conference on Communication, Control, and Computing, pages 1–9, 2004. [103] C. Radhakrishna Rao and Sujit Kumar Mitra. Generalized inverse of a matrix and its applications. In Proc. sixth Berkeley Symposium on Mathematical Statistics and Probability, pages 601–620, 1972. [104] M. G. Safonov and R. Y. Chiang. A Schur method for balanced- truncation model reduction. IEEE Transactions on Automatic Con- trol, 34(7), 1989. [105] M. G. Safonov, A. J. Laub, and G. Hartmann. Feedback properties of multivariable systems: The role and use of return difference matrix. IEEE Transactions on Automatic Control, 26(1):47–65, 1981. [106] P. Saint-Pierre. Approximation of the viability kernel. Applied Math- ematics and Optimization, 29(2):187–209, Mar 1994. 157 Bibliography [107] K.-H. Shim and M. E. Sawan. Singularly perturbed unified time sys- tems with low sensitivity to model reduction using delta operators. International Journal of Systems Science, 37(4):243–251, 2006. [108] Norihiko Shishido and Claire J. Tomlin. Ellipsoidal approximations of reachable sets for linear games. In Proc. IEEE Conference on Decision and Control, pages 999–1004, Sydney, Australia, December 2000. [109] O. Simanski, A. Schubert, R. Kaehler, M. Janda, J. Bajorat, R. Hof- mockel, and B. Lampe. Automatic drug delivery in anesthesia: From the beginning until now. In Proc. Mediterranean Conf. Contr. Au- tomation, Athens, Greece, 2007. [110] J. M. Siret, G. Michailesco, and P. Bertrand. Optimal approxima- tion of high-order systems subject to polynomial inputs. International Journal of Control, 26(6):963–971, 1977. [111] S. Skogestad and I. Postlethwaite. Multivariable Feedback Control; Analysis and Design. John Wiley & Sons, West Sussex, UK, 2007. [112] D. R. Smith. Decoupling and order reduction via the Riccati trans- formation. SIAM Review, 29(1):91–113, 1987. [113] D. M. Stipanović, I. Hwang, and C. J. Tomlin. Computation of an over- approximation of the backward reachable set using subsystem level set functions. In Proc. IEE European Control Conference, Cambridge, UK, September 2003. [114] G. Strang. Linear Algebra and Its Applications. Brooks Cole, 1988. [115] S. Syafiie, J. Niño, C. Ionescu, and R. De Keyser. NMPC for propofol drug dosing during anesthesia induction. In Nonlinear Model Predic- tive Control, volume 384, pages 501–509. Springer Berlin Heidelberg, 2009. [116] C. J. Tomlin, J. Lygeros, and S. Sastry. A game theoretic approach to controller design for hybrid systems. Proceedings of the IEEE, 88(7):949–970, 2000. [117] C. J. Tomlin, I. M. Mitchell, A. M. Bayen, and M. Oishi. Compu- tational techniques for the verification and control of hybrid systems. Proceedings of the IEEE, 91(7):986–1001, 2003. 158 [118] P. Varaiya. Reach set computation using optimal control. In Proc. KIT Workshop on Verification of Hybrid Systems. Verimag, Grenoble, 1998. [119] R. Vasudevan, V. Shia, Y. Gao, R. Cervera-Navarro, R. Bajcsy, and F. Borrelli. Safe semi-autonomous control with enhanced driver mod- eling. In Proc. American Control Conference, Montreal, QC, 2012. (to apprear). [120] R. Vidal, S. Schaffert, J. Lygeros, and S. Sastry. Controlled invariance of discrete time systems. In B. Krogh and N. Lynch, editors, Hybrid Systems and Control, LNCS 1790, pages 437–451, Berlin Heidelberg, 2000. Springer-Verlag. [121] H. Yazarel and G. J. Pappas. Geometric programming relaxations for linear system reachability. In Proc. American Control Conference, pages 553–559, Boston, MA, June 2004. [122] C. C. Zervos and G. A. Dumont. Multivariable self-tuning control based on Laguerre series representation. Adaptive Control Strategies for Industrial Use, 137:44–57, 1989. [123] K. Zhou, J. C. Doyle, and K. Glover. Robust and Optimal Control. Prentice Hall, Englewood Cliffs, NJ, 1996. 159 Appendix A Supplementary Materials for Schur-Based Decomposition A.1 On Assumption 3.1 (and 4.2) Proposition A.1. With A ∈ Rn×p and B ∈ Rm×p, the equation XA = B has a solution w.r.t. X ∈ Rm×n if and only if C (BT) ⊆ C (AT). Proof. By taking the transpose of both sides we have ATXT = BT. Let Y := XT. The equation ATY = BT is known to have at least one solution w.r.t. Y if and only if the column space of BT is in the image (or range) of the linear operator AT, i.e. C (BT) ⊆ C (AT). To check if the condition holds, we can check if span{~b1, . . . ,~bm} ⊆ span{~a1, . . . ,~an}, (A.1) where ~bi and ~ai denote the ith row of the matrices B and A, respectively. Or equivalently, we can check if the following rank condition holds: rank ([ AT|BT]) = rank (AT) . (A.2) A.2 Proof of Proposition 3.4 To prove Proposition 3.4 let us first state a simple lemma. Lemma A.1. Consider the backward reachable tube Reach[[0,t](K,U) of (2.12) over the interval [0, t], t ∈ [0, τ ], for a fixed τ ∈ R+. Denote by Reach[A[0,t](K) 160 A.2. Proof of Proposition 3.4 the backward reachable tube of its corresponding autonomous system ẋ = Ax. (A.3) The following inclusions hold: K ⊆ Reach[[0,t](K,U) ⊆ Reach[A[0,t](K) ∀t ∈ [0, τ ]. (A.4) Proof. Assume, without loss of generality, that U is a compact hyper-rectangular subset of Rp such that U = ∏pi=1 Ui, ui ∈ Ui = [U i,U i], 0 ∈ Ui. Notice that the trajectories of the autonomous system (A.3) approach those of the con- trolled system (2.12) as ξ := sup{‖u‖ : u ∈ U} (A.5) tends to zero. As such, we draw on the level set formulation of the backward reachable tube of system (2.12) and treat (A.3) as a particular form of (2.12) in which the control input u is diminished. It is well-known [89] that if K is represented as the zero sub-level set of some bounded and Lipschitz continuous implicit surface function g : Rn → R, i.e. K = {x | g(x) ≤ 0}, then the backward reachable tubeReach[[0,t](K,U) can be obtained as the zero sub-level set of the viscosity solution φ : Rn × [0, τ ]→ R of the modified terminal value HJB PDE ∇tφ(x, t) = −min {0, H (x,∇xφ(x, t))} , φ(x, τ) = g(x) (A.6) H(x, `) = sup u∈U 〈`, Ax+Bu〉 (A.7) with the Hamiltonian H(·, ·) and the costate vector `. Here, 〈·, ·〉 denotes the inner product. Thus, Reach[[0,t](K,U) = {x | φ(x, τ − t) ≤ 0}. The optimal Hamiltonian, in this case, can be determined analytically as H∗(x, `) = `TAx+ `TBu∗, u∗ = [u∗1 · · ·u∗p]T (A.8) 161 A.2. Proof of Proposition 3.4 with u∗i = U i if `Tbi < 0; [U i,U i] if `Tbi = 0; U i if `Tbi > 0 , i = 1, . . . , p (A.9) where bi is the i-th column vector of matrix B. Notice that the second term on the right hand side of (A.8) is always non-negative, i.e. `TBu∗ ≥ 0. Therefore we have ∇tφ(x, t) = 0 if `TAx ≥ −`TBu∗;|`TAx| − `TBu∗ otherwise. (A.10) When ξ = 0, the controlled system (2.12) is equivalent to the autonomous system (A.3) and the Hamiltonian (A.7) becomes H(x, `) = H∗(x, `) = `TAx. Consequently, (A.10) reduces to ∇tφ(x, t) ∣∣∣∣ ξ=0 =: ∇tφA(x, t) = 0 if `TAx ≥ 0;|`TAx| otherwise (A.11) where φA(·, ·) is to denote the implicit surface function whose zero sub-level set determines the backward reachable tube Reach[A[0,t](K) of (A.3). That is, Reach[A[0,t](K) = {x | φA(x, τ − t) ≤ 0}. Comparing (A.10) and (A.11) one can observe that not only the interval over which ∇tφ (the rate of surface change in time) is zero is shortened (i.e. `TAx ≥ 0 as opposed to `TAx ≥ −`TBu∗), but also its maximum (positive) value is increased (i.e. |`TAx| as opposed to |`TAx| − `TBu∗). Therefore, for all (x, t) ∈ Rn × [0, τ ] we have ∇tφA(x, τ − t) ≥ ∇tφ(x, τ − t) (A.12) =⇒ φA(x, τ − t) ≤ φ(x, τ − t) ≤ φ(x, τ) ≤ 0 (A.13) ⇐⇒ Reach[A[0,t](K) ⊇ Reach[[0,t](K,U) ⊇ K. (A.14) 162 A.2. Proof of Proposition 3.4 Notice that this result agrees with the intuitive interpretation that larger control authority (i.e. ξ 6= 0) implies a smaller unsafe minimal reachable tube. We are now ready to prove Proposition 3.4. Proof of Proposition 3.4. Using Lemma A.1 we have Reach[A[0,t](K) ⊇ Reach[[0,t](K,U) ∀t (A.15) where Reach[A[0,t](K) denotes the backward reachable tube of the autonomous system (A.3). Therefore, to prove Reach[[0,t](K,U) = K, ∀t ∈ [0, τ ], it is sufficient to show that Reach[A[0,t](K) = K, ∀t ∈ [0, τ ]. Let SΛS−1 be the eigen-decomposition of A. Conditions (ii) and (iii) imply Λ = λIn, λ ≥ 0. Rewriting the Hamiltonian of the HJB PDE (A.6) for the autonomous system (A.3) and using condition (i) we have H(x,∇xφA(x, t)) = 〈∇xφA(x, t), Ax〉 (A.16) = 〈∇xφA(x, t), SΛS−1x〉 (A.17) = λ〈∇xφA(x, t), x〉 ≥ 0 ∀(x, t) ∈ Rn × [0, τ ]. (A.18) The non-negativity of the Hamiltonian is due to the fact that K is convex and 0 ∈ K. Thus, the costate vector ∇xφA(x, t) at every point on the boundary constitutes an acute (hyper-) angle with respect to the trajectory x initiating from that point in forward time. This is schematically illustrated for a trivial planar system in Figure A.1. As a result, for all (x, t) ∈ Rn× [0, τ ] we have H(x,∇xφA(x, t)) ≥ 0⇐⇒ ∇tφA(x, t) = 0 (A.19) ⇐⇒ Reach[A[0,t](K) = K. (A.20) This concludes the proof. 163 A.3. Decomposed System Matrices for Examples 3.6.2 and 3.6.3 2A 2x 1x 1A 3A 3x K Figure A.1: Three sample costate vectors (`i) and trajectories (xi) initiating from the boundary of an arbitrarily-shaped convex target set K in the phase-plane of a simple planar system in forward time. Notice the non-negativity of 〈`i, xi〉 as shown in forward time. In backward time the trajectories are reversed and the eigenvalues are negated, hence the Hamiltonian is still non-negative. A.3 Decomposed System Matrices for Examples 3.6.2 and 3.6.3 A.3.1 4D Aircraft Dynamics (Example 3.6.2) The decomposed system matrices are: Ad = 0.1527 −0.1511 0.0312 0 0.1853 −0.1536 −0.0247 0.0065 0 0 −0.3169 −7.5973 0 0 0.1028 −0.4331 , Bd = 0 0 −0.1785 1.1598 . 164 A.3. Decomposed System Matrices for Examples 3.6.2 and 3.6.3 A.3.2 8D Distillation Column (Example 3.6.3) The system matrices after the first decomposition are: Ad = −0.6460 −2.7152 0.9186 −1.0340 −1.5499 0.0128 0.3436 0.0404 3.0705 −0.6461 −0.8642 0.6817 1.5878 −0.0215 −0.0567 −0.0336 0 0 −0.8627 1.6880 2.0531 0 0.0135 0 0 0 −0.6716 −0.8627 −0.5888 0 −0.2143 1.2635 0 0 0 0 −0.7357 −0.2275 −0.0082 −0.0021 0 0 0 0 0 −0.2259 0.0021 −0.0457 0 0 0 0 0 0 −0.0052 0.0024 0 0 0 0 0 0 0 −0.0755 Bd = 0 0 0 0 0 0 0 0 1.2886 −0.0504 0.3132 −0.2249 0.7117 −0.6994 0.0599 −0.3014 . The unidirectional coupling term in Ad is treated as a disturbance to the upper subsystem. The second level decomposition applied to each 4D subsystem results in Ad1 = −1.3594 −1.2819 0 0 1.0769 −0.3660 0 0 0 0 −0.9978 −2.8164 0 0 3.0041 −0.2943 , Bd1 = 0 0 0 0 0 0 0 0 , Ad2 = −0.0052 0.0026 0.0215 0.5192 0 −0.0755 −0.1334 0.1262 0 0 −0.7479 −0.1877 0 0 0.0339 −0.2137 , Bd2 = 0 0 0 0 1.3130 −0.0574 0.2260 −0.2934 . 165 Appendix B Supplementary Materials for Riccati-Based Decomposition B.1 Proofs of Propositions 4.2 and 4.3 Proof of Proposition 4.2. From the matrix inversion lemma, (Y+UCV )−1 = Y −1−Y −1U(C−1 +V Y −1U)−1V Y −1, with Y = −(δ+ 1)I, U = B1, C = I, and V = B†1 we have( B1B † 1 − (δ + 1)I )−1 = − 1 δ + 1 ( I + 1 δ B1B † 1 ) . (B.1) Using this, (4.14), (4.21), (4.23), (4.26), multiplicative and triangular in- equalities, and the fact that ‖B1B†1‖ ≥ 1 we obtain ‖δF (Z(δ))‖ ≤ |δ|(α(‖Z0‖+ ‖D‖)2 + β(‖Z0‖+ ‖D‖)) ≤ |δ| ( α ( ‖Z0‖+ 2‖A0‖‖Z0‖‖A0‖+ α‖Z0‖ )2 + β ( ‖Z0‖+ 2‖A0‖‖Z0‖‖A0‖+ α‖Z0‖ )) ≤ |δ|(9α‖Z0‖2 + 3β‖Z0‖) ≤ |δ| ( 9αγ2 ∥∥(B1B†1 − (δ + 1)I)−1∥∥2 + 3βγ ∥∥(B1B†1 − (δ + 1)I)−1∥∥) ≤ |δ| ( 9αγ2 ∣∣∣ 1 δ + 1 ∣∣∣2(1 + ∣∣∣1 δ ∣∣∣)2‖B1B†1‖2 + 3βγ ∣∣∣ 1 δ + 1 ∣∣∣(1 + ∣∣∣1 δ ∣∣∣)‖B1B†1‖) 166 B.2. Condition Number of the Modified Riccati Transformation = 1 |δ| ( |δ|+ 1 |δ + 1| )2 a+ ( |δ|+ 1 |δ + 1| ) b ∀δ ∈ R\{−1, 0}. Proof of Proposition 4.3. Notice from (4.25) and (4.28) that for large values of δ, Z can be closely approximated by its initial value Z0. Using (B.1), lim δ→±∞ ‖δF (Z(δ))‖ = lim δ→±∞ ∥∥∥ δ (δ + 1)2 Q1(I + 1 δ B1B † 1)P1Q1 × (I + 1 δ B1B † 1) + δ δ + 1 P2Q1(I + 1 δ B1B † 1) ∥∥∥ (B.2) = ‖0 + P2Q1‖ = ‖Γ‖ (B.3) with Q1 := ( B2B † 1A12 − A22 )−1 Γ, P1 := (A12 − B1B†1A12), and P2 := (B2B † 1A12 −A22). B.2 Formulating an Upper-Bound on the Condition Number of the Modified Riccati Transformation Lemma B.1. The condition number κ(T ) of the Riccati transformation matrix T = T1T2 is bounded by κ(T ) ≤ max{1 + µ, 1 + ν(1 + µ)} ·max{1 + ν, 1 + µ(1 + ν)} (B.4) with constants µ := 2‖N0‖‖M0‖ ‖N0‖+ ‖δF (Z)‖‖M0‖ + ‖M0‖, ν := ‖B2B†1‖+ (1 + ‖B1B†1‖) [ 2‖A0‖‖Z0‖ ‖A0‖+ ‖B1B†1A12 −A12‖‖Z0‖ + ‖Z0‖ ] . Proof. Let A = [Aij ] be any partitioned matrix. Then ‖A‖ ≤ ∥∥[‖Aij‖]∥∥. 167 B.3. Decomposed System Matrices for Example 4.5.3 and Section 4.5.5 Since T = [ I M −L I−LM ] and T−1 = [ I−ML −M L I ] , we have κ(T ) = ‖T‖‖T−1‖ ≤ ∥∥∥[ 1 ‖M‖‖L‖ 1+‖L‖‖M‖ ]∥∥∥ ∥∥∥[ 1+‖M‖‖L‖ ‖M‖‖L‖ 1 ]∥∥∥ . (B.5) From (4.10) we find ‖L‖ ≤ ‖B2B†1‖ + ‖Z‖ ( 1 + ‖B1B†1‖ ) . But from (4.23) and (4.26) we know that ‖Z‖ is bounded above by 2‖A0‖‖Z0‖‖A0‖+‖B1B†1A12−A12‖‖Z0‖+‖Z0‖, and similarly from (4.33) and (4.36), ‖M‖ is bounded above by 2‖N0‖‖M0‖ ‖N0‖+‖δF (Z)‖‖M0‖ + ‖M0‖. Therefore, ‖L‖ ≤ ν and ‖M‖ ≤ µ and con- sequently, κ(T ) ≤ max{1 + µ, 1 + ν(1 + µ)} ·max{1 + ν, 1 + µ(1 + ν)}. B.3 Decomposed System Matrices for Example 4.5.3 and Section 4.5.5 B.3.1 Arbitrary 6D System (Example 4.5.3) The decomposed system matrices are: A′′ = 3.3126 0.7676 2.4511 0 0 0 0.6223 −1.5072 −2.3105 0 0 0 −0.0989 −0.3285 0.6852 0 0 0 0.0944 0.0042 −0.1299 0.1880 0.0445 0.3528 −0.0540 −0.4170 −0.0461 0.2802 0.0888 0.1593 −0.1474 0.1949 0.2443 0.0848 0.0888 0.1197 , B′′ = 0.1469 0.2657 −0.7988 2.4582 −2.3854 −0.3955 0 0 0 0 0 0 . 168 B.3. Decomposed System Matrices for Example 4.5.3 and Section 4.5.5 B.3.2 Maximal Reachability Example (Section 4.5.5) The decomposed system matrices are: A′′ = 0.5912 −0.2477 0 0 −0.4583 −0.2017 0 0 0.0197 0.1351 −0.1905 −0.1878 −0.0825 −0.1479 0.0012 −0.0633 , B′′ = 0.6846 2.4813 0 0 . 169 Appendix C Other Backward Reachability Constructs1 Some additional backward constructs formed under the constrained dynam- ical system (2.1), their connections to one another and to the constructs used within this thesis (see Chapter 2) are presented here, aiming to help the reader attain a more complete picture. C.1 Definitions and Connections Definition C.1 (Minimal Reachable Set). The minimal reachable set at time t is the set of initial states such that, for every input u(·), the trajecto- ries emanating from those states reach K exactly at time t: Reach[t(K,U) := { x0 ∈ X | ∀u(·) ∈ U[0,t], xux0(t) ∈ K } . (C.1) Definition C.2 (Invariance Kernel). The (finite-horizon) invariance kernel of K is the set of initial states in K such that the trajectories emanating from those states remain within K for all time t ∈ [0, τ ] for all input u(·): Inv[0,τ ](K,U) := { x0 ∈ X | ∀u(·) ∈ U[0,τ ], ∀t ∈ [0, τ ], xux0(t) ∈ K } . (C.2) Definition C.3 (Continual Reachable Set [59, 60]). The continual reach- able set defined over the time horizon [0, τ ] is the set of initial states in K for which, for any given time t ∈ [0, τ ], there exists a u(·) such that the 1A part of this chapter is based on [59, 60]. The main results of these papers, however, have not been presented in this thesis. 170 C.1. Definitions and Connections trajectories emanating from those states reach K at t: Reachγ[0,τ ](K,U) := { x0 ∈ X | ∀t ∈ [0, τ ], ∃u(·) ∈ U[0,t], xux0(t) ∈ K } . (C.3) The following inclusions are complementary to V iab[0,τ ](K,U) ⊆ K ⊆ Reach[[0,τ ](K,U) ⊆ Reach][0,τ ](K,U) described in Proposition 2.1. Proposition C.1. Inv[0,τ ](K,U) ⊆ V iab[0,τ ](K,U) ⊆ Reachγ[0,τ ](K,U) ⊆ K. (C.4) Proof. That Inv[0,τ ](K,U) ⊆ V iab[0,τ ](K,U) is well-known [5]. To show V iab[0,τ ](K,U) ⊆ Reachγ[0,τ ](K,U), take x0 ∈ V iab[0,τ ](K,U). Therefore, ∃u(·) ∈ U[0,τ ] ∀t ∈ [0, τ ] xux0(t) ∈ K =⇒ ∀t ∈ [0, τ ] ∃u(·) ∈ U[0,t] xux0(t) ∈ K ⇐⇒ x0 ∈ Reachγ[0,τ ](K,U). To show Reachγ[0,τ ](K,U) ⊆ K, take x0 ∈ Reachγ[0,τ ](K,U) and let τ = 0. x0 must also belong to K. The maximal reachable tube and the invariance kernel are duals of one another as mentioned in Section 1.2: Proposition C.2 ([18], [79]). Reach][0,τ ](Kc,U) = (Inv[0,τ ](K,U))c. (C.5) Unlike the maximal reachable tube and sets, the minimal reachable tube cannot be constructed from the union of the minimal reachable sets as men- tioned in Section 2.2: Proposition C.3 ([86]). Reach[[0,τ ](K,U) ⊇ ⋃ t∈[0,τ ] Reach[t(K,U). (C.6) Among Lagrangian methods, the technique in [72] has been extended to handle universally quantified inputs. Therefore, it is also capable of computing the minimal reachable sets. As a by-product of this feature, the same technique can also be used to directly compute the invariance kernel. 171 C.1. Definitions and Connections Proposition C.4 ([59, 79]). Inv[0,τ ](K,U) = ⋂ t∈[0,τ ] Reach[t(K,U). (C.7) Proof. x0 ∈ ⋂ t∈[0,τ ]Reach [ t(K,U) ⇐⇒ ∀t ∈ [0, τ ] ∀u(·) ∈ U[0,τ ] xux0(t) ∈ K ⇐⇒ ∀u(·) ∈ U[0,τ ] ∀t ∈ [0, τ ] xux0(t) ∈ K ⇐⇒ x0 ∈ Inv[0,τ ](K,U). This can also be verified from (C.5) and (2.9) and the simple fact that Reach[t(K,U) = (Reach]t(Kc,U))c. Finally, the continual reachable set can be expressed in terms of the maximal reachable sets as: Proposition C.5 ([59]). Reachγ[0,τ ](K,U) = ⋂ t∈[0,τ ] Reach]t(K,U). (C.8) Proof. x0 ∈ Reachγ[0,τ ](K,U) ⇐⇒ ∀t ∈ [0, τ ] ∃u(·) ∈ U[0,t] xux0(t) ∈ K ⇐⇒ ∀t ∈ [0, τ ] x0 ∈ Reach]t(K,U)⇐⇒ x0 ∈ ⋂ t∈[0,τ ]Reach ] t(K,U). 172
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Scalable techniques for the computation of viable and...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Scalable techniques for the computation of viable and reachable sets : safety guarantees for high-dimensional… Kaynama, Shahab 2012
pdf
Page Metadata
Item Metadata
Title | Scalable techniques for the computation of viable and reachable sets : safety guarantees for high-dimensional linear time-invariant systems |
Creator |
Kaynama, Shahab |
Publisher | University of British Columbia |
Date Issued | 2012 |
Description | Reachability analysis and viability theory are key in providing guarantees of safety and proving the existence of safety-preserving controllers for constrained dynamical systems. The minimal reachable tube and (by duality) the viability kernel are the only constructs that can be used for this purpose. Unfortunately, current numerical schemes that compute these constructs suffer from a complexity that is exponential in the dimension of the state, rendering them impractical for systems of dimension greater than three or four. In this thesis we propose two separate approaches that improve the scalability of the computation of the minimal reachable tube and the viability kernel for high-dimensional systems. The first approach is based on structure decomposition and aims to facilitate the use of computationally intensive yet versatile and powerful tools for higher-dimensional linear time-invariant (LTI) systems. Within the structure decomposition framework we present two techniques – Schur-based and Riccati-based decompositions – that impose an appropriate structure on the system which is then exploited for the computation of our desired constructs in lower-dimensional subspaces. The second approach is based on set-theoretic methods and draws a new connection between the viability kernel and maximal reachable sets. Existing tools that compute the maximal reachable sets are efficient and scalable with polynomial complexity in time and space. As such, these scalable techniques can now be used to compute our desired constructs and therefore provide guarantees of safety for high-dimensional systems. Based on this new connection between the viability kernel and maximal reachable sets we propose a scalable algorithm using ellipsoidal techniques for reachability. We show that this algorithm can efficiently compute a conservative under-approximation of the viability kernel (or the discriminating kernel when uncertainties are present) for LTI systems. We then propose a permissive state-feedback control strategy that is capable of preserving safety despite bounded input authority and possibly unknown disturbances or model uncertainties for high-dimensional systems. We demonstrate the results of both of our approaches on a number of practical examples including a problem of safety in control of anesthesia and a problem of aerodynamic flight envelope protection. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2012-08-01 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0072951 |
URI | http://hdl.handle.net/2429/42856 |
Degree |
Doctor of Philosophy - PhD |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2012-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2012_fall_kaynama_shahab.pdf [ 3.58MB ]
- Metadata
- JSON: 24-1.0072951.json
- JSON-LD: 24-1.0072951-ld.json
- RDF/XML (Pretty): 24-1.0072951-rdf.xml
- RDF/JSON: 24-1.0072951-rdf.json
- Turtle: 24-1.0072951-turtle.txt
- N-Triples: 24-1.0072951-rdf-ntriples.txt
- Original Record: 24-1.0072951-source.json
- Full Text
- 24-1.0072951-fulltext.txt
- Citation
- 24-1.0072951.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-0072951/manifest