Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Scalable techniques for the computation of viable and reachable sets : safety guarantees for high-dimensional.. Kaynama, Shahab 2012-12-31

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

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

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 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 underii  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 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.  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 Schurbased 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 reachability analysis of linear time-invariant systems,” in Proc. IEEE Conference 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. Dumont, “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 appeared in: iv  Preface [v] S. Kaynama, M. Oishi, I. M. Mitchell, and G. A. Dumont, “The continual reachability set and its computation using maximal reachability techniques,” in Proc. Joint IEEE Conference on Decision and Control, 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, “Fixedcomplexity piecewise ellipsoidal representation of the continual reachability 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 simulations. I was responsible for the bulk of the writing of the papers. My supervisors 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 continuoustime 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 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  ii  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  iv  Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . .  vi  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  x  List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  xi  Abstract Preface  List of Tables  . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii  Acknowledgments 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  . . . . . . . . . . . . . . . . . . . . . . . . . .  8  1.5  Related Work 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 . . . . . . . . . . .  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  Approach II: Set-Theoretic Methods . . . . . . . . . . . . . .  25  2.5.1  Objective and Problem Statement . . . . . . . . . . .  25  2.5.2  Definitions and Preliminaries . . . . . . . . . . . . . .  25  2 Preliminaries and Problem Formulation  2.5  2.6  Robust Reachability Analysis: Competing Inputs  . . . . . .  26  2.7  To Over- or Under-Approximate?  . . . . . . . . . . . . . . .  29  . . . . . . . . . . . .  30  3.1  Disjoint Control Input . . . . . . . . . . . . . . . . . . . . . .  31  3.2  Non-Disjoint Control Input . . . . . . . . . . . . . . . . . . .  33  3.3  Reachability in Lower Dimensions  34  3 Schur-Based Structure Decomposition  3.3.1  . . . . . . . . . . . . . . .  Formulating an Upper-Bound on Growth of  1 Z[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  Summary and Further Discussions . . . . . . . . . . . . . . .  48  4 Riccati-Based Structure Decomposition . . . . . . . . . . . .  51  3.7  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  . . . . . . . . . . . . . . . . . . . .  4.4  Reachability in Lower Dimensions  . . . . . . . . . . . . . . .  65  . . . . . . .  67  4.4.3  Conservatism Due to Projection . . . . . . . . . . . .  67  4.4.4  Implications of Riccati-Based Reachable Tube  . . . .  68  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  Summary and Further Discussions . . . . . . . . . . . . . . .  79  5 Set-Theoretic Methods: Lagrangian Tools for Viability . .  81  4.6  5.1  5.2  Formulating an Upper-Bound on Growth of  4.4.2  Dimension vs. Magnitude of Uncertainty  63  . .  4.5  4.4.1  2 Z[0,τ ]  62  The Viability Kernel in Terms of Maximal Reach Sets . . . .  82  5.1.1  Continuous-Time Systems  . . . . . . . . . . . . . . .  82  5.1.2  Discrete-Time Systems  . . . . . . . . . . . . . . . . .  87  . . . . . . . . . . . . . . . . . . .  88  Computational Algorithms 5.2.1  5.3  5.4  A Piecewise Ellipsoidal Approach  . . . . . . . . . . .  89  . . . . . . . . . . . . . . . . . . . . . . .  97  5.3.1  Flight Envelope Protection (Continuous-Time) . . . .  97  5.3.2  Safety in Anesthesia Automation (Discrete-Time)  Practical Examples  Summary and Further Discussions . . . . . . . . . . . . . . . 103  6 Robustness and Safety-Preserving Control Synthesis 6.1  Discriminating Kernel vs. Robust Maximal Reach Sets 6.1.1  6.2  6.3  . . 100  . . . 106 . . . 106  Under-Approximating the Discriminating Kernel . . . 107  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  Numerical Examples . . . . . . . . . . . . . . . . . . . . . . . 122  viii  Table of Contents  6.4  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  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) A.3.2 8D Distillation Column (Example 3.6.3) B Supplementary Materials for Chapter 4 B.1 Proofs of Propositions 4.2 and 4.3  . . . . . . . . 164 . . . . . . . 165  . . . . . . . . . . . 166  . . . . . . . . . . . . . . . 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 C.1 Definitions and Connections  . . . . . . . . . . 170  . . . . . . . . . . . . . . . . . . 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 Reacht (K, U) . . . . . . . . . . . .  16  The maximal reachable tube Reach[0,τ ] (K, U) . . . . . . . . .  17  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  2.2 2.3 2.4  The minimal reachable tube Reach[0,τ ] (K, U) . . . . . . . . .  18  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 systems 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. . . . . . . . . . . .  5.3  A finer time interval partition results in better approximation 95  5.4  Convergence plot of the error as a function of |P | . . . . . . .  5.5 5.6 5.7  Under-approximation of V iab[0,2] (K, U) for Example 5.3.1 . .  Under-approximation of V iab[0,2] (K, U) for Example 5.3.1 . .  86 95 98 99  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 number 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ˇsan Stipanovi´c (U. of Illinois 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 interxiii  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 education 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´e 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 #CHRPJ-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 facilitating 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 infusion 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 enforce 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 thorough treatment of such techniques as well as other application areas where 1  1.4 Response  1 Setpoint  0.6  Safety bounds  0.2 0  1000  2000  3000 Time (s)  4000  5000  Input (Infusion Rate)  Output (Pseudo-occupancy)  1.1. Motivation  0.12  0.08  0.04  0  (a) Response of the closed-loop system.  1000  2000  3000 Time (s)  4000  5000  (b) Synthesized control law.  Figure 1.1: Control of drug plasma concentration in a high-responder patient (#80) undergoing 1.5 hr surgery using an inadequatelydesigned MPC. The input (drug infusion rate) is bounded above and below (u ∈ [0, 0.8]). Clearly, the hard constraint on the output (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 underwater 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 control 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 uncertain oscillatory systems [23] and verification of temporal properties of a toggle circuit [44], etc. Reachability analysis identifies the set of states forward (backward) reachable 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 trajectories 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. Similarly, 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 introduced in [86]. Their corresponding constructs differ in how the time variable 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 viability 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 reachability 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 impossible. 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 constraints.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 ability to compute the minimal reachable tube and the viability kernel can be paramount: Guarantees of safety and synthesis of safety-preserving controllers 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, 1 The 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. 2 One 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 computation of the controlled-invariant subset (viability kernel) when all constraints are polytopes. However, similarly to Eulerian methods, this technique suffers from an exponential complexity 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 higherdimensional 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 nonconvex 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 discretetime (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 techniques 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 conditions) 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 lowerdimensional 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 proposed techniques are of primary benefit to the computation of the minimal reachable tube (or by duality, the viability kernel) for which Eulerian methods 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 safetypreserving control laws.  7  1.5. Related Work  1.5  Related Work  Complexity reduction for reachability analysis has been addressed by a number of researchers. Methods to compute the reachability constructs (maximal or minimal) for higher-dimensional systems can be divided into three main categories. First are techniques that take advantage of certain representations of sets [30, 40, 64, 65, 74, 76, 78, 108]. Second are techniques that make use of model approximation [42, 47], hybridization [3], projection [44, 91] and structure decomposition [48, 113, 121]. Finally, third are methods that combine the approaches from the first two categories; For example, 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 reduction 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 differential equations (PDEs) over-approximates the projection of the actual minimal reachable tube in lower dimensional subspaces, with the unmodeled dimensions 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 viscosity 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 control set is a two-dimensional, axis-aligned rectangle), the minimal reachable tube equals the intersection of two reachable tubes each with a constant input. 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 x˙ = 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 Lf B(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 Lf B(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 computed set (viability kernel or minimal reachable tube) by choosing the safetypreserving 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 viability 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 satisfaction of the closed-loop system and therefore the generated control laws can be regarded as safety-preserving. In both cases discussed above the computational complexity of synthesizing such control laws is at best equal to that of computing the corresponding viability kernel or minimal reachable 3  This 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). 4 While 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 certificate technique. The technique has been extended in [102] to synthesize safety-preserving controllers using density functions and convex optimization methods for polynomial nonlinear systems with semi-algebraic constraints 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 behavior 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 bisimilar 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 thesis. 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 reachability analysis. These two chapters are based on the assumption that the system dynamics are LTI and aim to facilitate the applicability of computationally intensive Eulerian methods for higher-dimensional LTI systems. Chapter 5 describes a set-theoretic approach that enables the use of efficient Lagrangian methods for the computation of the viability kernel, circumventing entirely the need to employ computationally intensive Eulerian methods. An efficient algorithm is presented that approximates the viability 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. Chapter 7 summarizes the thesis reiterating its major contributions, and provides directions for future research. Finally, supplementary materials are provided in the Appendices: Appendix 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=0  Av = max 1≤j≤n v  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 inf d(x, y), sup inf d(x, y) . x∈X y∈Y  (1.3)  y∈Y x∈X  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 cl X := 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 Reacht (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):  Reacht (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) 1 Also 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 2 Also 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). 3 The 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 properties 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 show K ⊆ 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 “safetypreserving control synthesis”. To see this let us begin with maximal reachability 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,τ ]  = t∈[0,τ ]  x0 ∈ X | ∃u(·) ∈ U[0,t] , xux0 (t) ∈ K Reacht (K, U).  (2.8) (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 interchanged. This connection, which has been recognized and extensively used since the earlier works in the literature (e.g. [20, 22, 69]), was also formalized 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 handle universally quantified inputs. However, since the time variable is quantified 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 safetypreserving 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, 4 Under-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 decomposition techniques that will be presented in Chapters 3 and 4, aims to enable the use of powerful Eulerian methods on higher-dimensional LTI systems. 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 2.4.1  Approach I: Structure Decomposition Objective and Problem Statement  Consider the case in which (2.1) is a continuous-time LTI system x˙ = 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 reachability 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 −1 AT T −1 B . A linear transformation of a set V ⊆ Rn under the same mapping is defined as Y = T −1 V := {y ∈ Rn | y = T −1 v, v ∈ V}.  22  2.4. Approach I: Structure Decomposition Definition 2.5 (Unidirectionally Coupled). The LTI system that consists of two subsystems x˙ 1 = A1 x1 + ∆ x2 ,  (2.14)  x˙ 2 = A2 x2  (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 characterized by the maximum row sum ∆ . Definition 2.6 (Unidirectionally Weakly-Coupled). Let there be a nonT T singular transformation matrix T ∈ Rn×n , such that [z1T , z2T ]T = T −1 [xT 1 , x2 ] ,  and  z˙1 = A1 z1 + ∆ z2  (2.16)  z˙2 = A2 z2 .  (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 −1 x in which (2.12) can be partitioned into N subsystems as  z˙i = Ai zi + gi (u),  i = 1, . . . , N,  (2.19)  with gi (u) := Bi u. The input u = [u1 , . . . , up ]T ∈ U ⊂ Rp is disjoint across these subsystems if ∀s ∈ {1, . . . , p}, ∀i, j ∈ {1, . . . , N }, i = j, ∂gj (u) ∂gi (u) =0 → = 0. ∂us ∂us  (2.20)  23  2.4. Approach I: Structure Decomposition Example 2.1. Consider the subsystem trajectories of z˙i = Ai zi + Bi u with Bi = [Bi1 , Bi2 , Bi3 ], u = [u1 , u2 , u3 ]T , and i ∈ {1, 2}. The input u is disjoint if for example z1 (t) z2 (t)  =  eA1 (t−t0 )  0  z1 (t0 )  0  eA2 (t−t0 )  z2 (t0 )  t  +  eA1 (t−r)  0  B11  0  eA2 (t−r)  0  t0  0    u1 (r)   u2 (r) dr. (2.21) u3 (r)  0  B22 B23  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)  =  eA1 (t−t0 )  0  z1 (t0 )  0  eA2 (t−t0 )  z2 (t0 )   t  + t0  eA1 (t−r)  0  0  eA2 (t−r)  B11 B12 0  0  B22 B23  u1 (r)       u2 (r)  dr. (2.22) u3 (r)  Definition 2.8 (ETUC). A subsystem is said to be externally-triviallyuncontrollable (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 (F T ⊗ Ik ) + (Im ⊗ E) − vec(H) = rank (F T ⊗ 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 ) = 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 U T M U = 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 = M11 M12 0 M22  . 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 2.5.1  Approach II: Set-Theoretic Methods Objective and Problem Statement  For the generic, continuously-valued system (2.1) we define our second objective as follows. Problem 2.2 (Objective II). Express the viability kernel V iab[0,τ ] (K, U) in terms of the maximal reachable sets Reacht (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 x˙ = 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 |, |P |  • the size of the largest interval by P := maxk=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 functions {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 concatenation 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 x − a . a∈A  (2.28)  For a fixed set A, the map x → 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 x(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,v x0 : [0, t] → X that satisfies the initial condition  xu,v x0 (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 specified. The control input follows a feedback strategy, i.e. u(t) = u ˆ(x(t), t) ∈ U  ∀t ∈ [0, τ ]. We assume non-anticipative strategies for the disturbance input. 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: Reacht (K, U, V) := {x0 ∈ X | ∃u(·) ∈ U[0,t] ,  ∀ρ : U[0,t] → V[0,t] , xu,ρ[u] (t) ∈ K}. (2.31) x0  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,τ ] ,  (t) ∈ K}. (2.32) ∀u(·) ∈ U[0,τ ] , ∃t ∈ [0, τ ], xu,ρ[u] x0  Definition 2.15 (Discriminating Kernel). The (finite-horizon) discriminating 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] (t) ∈ K . (2.33) x0  Note that with V = {0} the discriminating kernel reduces to the viability  kernel under the deterministic system x˙ = 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 5  The 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.7  (2.35)  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 underapproximation 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 underapproximated 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 viability 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 inspired by a model reduction algorithm for systems with unstable modes [82, 110]. Our method decomposes LTI systems into either completely decoupled or weakly-coupled subsystems. Reachability analysis can be performed on each resulting subsystem independently. Back projecting and intersecting each of the lower-dimensional reachable tubes provides an overapproximation 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 underapproximation 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 transformation 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 decom1  A 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 isolated subsystem, whereas for the remaining subsystem, the effect of coupling is accounted for by treating the coupling terms as disturbance and performing reachability with competing inputs. For both decoupled and unidirectionally weakly-coupled decompositions, the intersection of back projections of the lower dimensional reachable tubes is an over-approximation of the actual reachable tube in the transformed coordinate space. When the control input across the decomposed subsystems is non-disjoint, a constrained optimization 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 decomposition) 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  A11 A12 B1 0  A22 B2  (3.1)  with A11 ∈ Rk×k , A12 ∈ Rk×(n−k) , A22 ∈ R(n−k)×(n−k) , B1 ∈ Rk×p , and  B2 ∈ 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  A11 X − X A22 + A12 = 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 A11 A11 X−X A22 +A12 B1 A22  0  B2  =  A11 0  0  B1  A22 B2  .  (3.4)  Notice that the resulting subsystems  A11 B1 and A22 B2 have been effectively decoupled through the coordinate transformation z = T −1 x, T = U W . As we shall see in Section 3.3, reachability analysis (in this transformed 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)  A11 Q − QA22 + A12  (3.5)  results in unidirectionally weakly-coupled subsystems w.r.t. (3.1). Proof. Consider Ac := A11 X − X A22 + A12 = 0 in (3.4). It is clear that in  the transformed coordinate space characterized by z = (U W )−1 x, the state vector z2 evolves independently of z1 since z˙2 = A22 z2 + B2 u2 . However, z1 is affected by z2 through Ac . That is, we have z˙1 = A11 z1 + B1 u1 + Ac z2 . (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{ A11 Q − QA22 + A12 Rk×(n−k) }. imply that  Then the hypothesis A12 < A11  X∗  X∗  −  X ∗A  22  | Q ∈  + A12 would  = 0 can never be a solution. Since there are no constraints in  (3.5) imposing this restriction, by contradiction we conclude that A11 X ∗ −  X ∗ A22 + A12 ≤ A12 . 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 infinity 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 (A22 , B2 ) or the pair (A11 , B1 ) 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 (B1T ) ⊆ C (B2T ) with C (·) the column-space operator. Proposition 3.3. The transformation (3.3) with A11 Q − QA22 + A12  X = arg min Q∈Rk×(n−k)  (3.6)  subject to QB2 = B1 results in unidirectionally coupled subsystems. Moreover, (A11 , B1 ) is ETUC. Proof. Assumption 3.1 is the necessary and sufficient condition for solvability of the overdetermined equality constraint in (3.6) (cf. Appendix A.1). To see the trivial-uncontrollability of (A11 , B1 ) consider B := W −1 B in (3.4). We have B1 B2  :=  Ik 0  −X  In−k  B1 B2  =  B1 − X B2 B2  .  (3.7)  Constraining the optimizer in (3.6) to choose from the class of solutions X ∈ Rk×(n−k) | X B2 = B1 simply enforces B1 = 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 coordinate space z = T −1 x, T := U W 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 −1 K 2:  Zτi ← ProjSi (Zτ ), ∀i ∈ {1, 2}  project onto ith subspace  For lower subsystem: 3:  2 2 Z[0,τ ] ← Reach[0,τ ] (Zτ , U)  For upper subsystem: 4:  Treat Ac z2 as disturbance  5:  ζ ← z2 ≡ supv∈Z 2  [0,τ ]  6: 7: 8:  v  Ac := A11 X − X A22 + A12  Compute upper-bound Ac z2 ≤ Ac ζ consrv.  1 1 Z[0,τ ] ←− Reach[0,τ ] (Zτ , U, B( Ac ζ))  1 , Z2 ) return(Z[0,τ ] [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: 1 2 S1 (Proposition 3.1 is used): For both Z[0,τ ] and Z[0,τ ] , the single input  is control and it is universally quantified.  1 S2 (Proposition 3.2 is used): For Z[0,τ ] the control input is universally  quantified while the disturbance input (unidirectional coupling) is ex-  2 istentially quantified. For Z[0,τ ] the single input is control and it is  universally quantified.  1 S3 (Proposition 3.3 is used): For Z[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 al2 gorithm is effectively replaced by Reach[0,τ ] (Zτ1 , B( Ac ζ)). For Z[0,τ ]  the single input is control and it is universally quantified.  The over-approximation of the actual minimal reachable tube of the fullorder system in X can be obtained using the following lemma.  35  3.3. Reachability in Lower Dimensions i Lemma 3.1 ([91, 113]). Let Z[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): 2 1 T (Z[0,τ ] × S2 ) ∩ (S1 × Z[0,τ ] ) ⊇ Reach[0,τ ] (K, U).  3.3.1  (3.9)  1 Formulating an Upper-Bound on Growth of Z[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 quantified. Consequently, this disturbance together with the dynamics strive to enlarge the reachable (unsafe) set as much as possible. This allows us to formulate 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 ∈ Zτ1 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−A11 τ z¯1 − 1 z1 ∈ Z[0,τ ].  τ  e−A11 (τ −r) ϑ(r)dr,  (3.10)  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−A11 τ z¯1 ≤ =  τ  e  A11 (τ −r)  e  A11 τ  A11  −1 M  = ≤  Ac ζdr  (3.12)  0  lim  M →∞  i=1 M  lim  M →∞  i=1  Ac ζ  (3.13)  τ i A11 i−1 Ac ζ i! √ i−1 τ i k σ(A11 ) Ac ζ i!  =: µ[0,τ ]  (3.14)  (3.15) (3.16)  where σ(·) is the largest singular value operator, and k is the dimension of 1 the ETUC subsystem. Therefore, an upper-bound for how much Z[0,τ ] can  grow in backward time can be written as  1  Z[0,τ ] ⊆   t∈[0,τ ]  e−A11 t Zτ1  ⊕ 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∈Z 2  [0,τ ]  v , and the largest  singular value of the upper subsystem σ(A11 ) can all affect the conservatism 1 . Moreover, given k and τ , the flexibility of the of the reachable tube Z[0,τ ]  Schur form in placing the eigenvalues in any order along the block-diagonals of A can be exploited to make this subsystem evolve with slower dynamics. Through various tests we were able to confirm that doing so could potentially 1 −A11 t and µ prevent the excessive growth of Z[0,τ [0,τ ] . ] by influencing both e  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 antistable 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  λi = λj K  nonconvex set  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 violate 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 fullorder 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 x˙ = 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 A− 0  Ac  0  A+ B2  (3.19)  partitioned such that A+ and A− contain only non-negative and strictlynegative eigenvalues, respectively. If A+ and K satisfy the conditions (i), 39  3.5. Extension to Switched Linear Systems  K  λi = λj Algebraic Geometric multiplicity > multiplicity  λi = / λj  Re(λ) = 0 ∧ Im(λ) = / 0  Im(λ) = / 0 ∧ Re(λ) = / 0  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. Reachability analysis is performed only for the upper subsystem resulting in further reduction of complexity by avoiding altogether the reachable tube computation 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 2 −1 K), U) ≡ be calculated according to (3.9) with Z[0,τ ] = Reach[0,τ ] (ProjS2 (T 1 −1 K), B( A ζ)). ProjS2 (T −1 K) and Z[0,τ c ] = Reach[0,τ ] (ProjS1 (T  Note that linear transformation preserves convexity. Therefore the pro-  jection of the transformation of K onto the lower subspace (i.e. Zτ2 :=  ProjS2 (T −1 K)) 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 dynamical 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) → Ai x + Bi u,  (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  Tj Reach[0,τ ] qj , Tj−1 Ti Reach[0,τ ] qi , Ti−1 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 accommodate both existentially and universally quantified inputs, we demonstrate 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 complexity 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   −0.5672 −0.7588 −0.6282   A =  3.1364 −1.1705 2.3247 , 1.8134 −1.7689 −2.6930     0.0731 −0.1639   B = −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 −1 K 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:   −1.6653 −3.4560  0      T −1 AT =  0   1.8706 −1.4653 , 0 0 −1.3000 Hence, the decoupled subsystems are z˙1 =    −0.7530  0      T −1 B =  0   0.0640 . 0 0.2500 −1.6653 −3.4560 1.8706 −1.4653  z1 +  −0.7530 0.0640  u1  and z˙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  (c) (z1 , z3 ) cross-section  (b) (z1 , z2 ) 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 projections) were 5823.73 s and 22.87 s, respectively.  43  3.6. Numerical Examples  3.6.2  4D Aircraft Dynamics  Consider longitudinal aircraft dynamics x˙ = Ax + Bδe ,   −0.0030 0.0390 0 −0.3220   −0.0650 −0.3190 7.7400 0  ,  A=  0.0200 −0.1010 −0.4290 0   0 0 1 0   0.0100   −0.1800   B=  −1.1600   0   ˙ θ]T ∈ R4 comprised of deviations in aircraft velocity with state x = [u, v, θ, [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 −1 x, x ∈ K} where T is the transformation matrix obtained through our method. We first decompose 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 Appendix 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 reduction. 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 2  crad = 0.01 rad ≈ 0.57◦ .  44  3.6. Numerical Examples  Figure 3.5: Schur-based (solid dark) vs. actual (transparent light) finitehorizon 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.   −0.5774  −2.7290   0   0  A=  0   0   0  0  3.0567 −0.7147 0 0 0 0 0 0  0.0073 −0.3430 −0.3891 1.3640 0 0 0 0  −0.8121 1.5321 −0.9956 −1.3363 0 0 0 0  0.3034 0.6643 0.0182 −0.9037 −0.7357 0 0 0  −0.3035 0.2896 0.0235 −0.4686 −0.2275 −0.2259 0 0  0.0072 −0.0013 0.0049 −0.0009 −0.0082 0.0021 −0.0052 0   −0.1542  0.0926  0.0506  −0.1887   −0.0021  −0.0457   0.0024 −0.0755  −0.0335 −0.1228  −0.4534 −0.0711  −0.8005 −0.2612  0.5497 −0.1344  1.2886 −0.0504  0.3132 −0.2249  0.7117 −0.6994  0.0599 −0.3014  B=  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 w1 = [w1 , w2 ]T , w2 = [w3 , w4 ]T , q1 = [q1 , q2 ]T , and q2 = [q3 , q4 ]T . Notice that R4 T2−1 z1 , and R8  q = [q1T , q2T ]T = T3−1 z2 , R4  w = [w1T , w2T ]T =  z = [z1T , z2T ]T = T1−1 x with z1 , z2 ∈ R4 .  We assume that the target (unsafe) set K ⊂ R8 is chosen such that  the transformations T1−1 ∈ R8×8 , T2−1 ∈ R4×4 , and T3−1 ∈ 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 of Wτ 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  20 w4  w2  20 0  −20  −20 −20  0 w1  20  −20  0 w3  20  −20  0 q3  20  20 q4  20 q2  0  0  −20  0  −20 −20  0 q1  20  Figure 3.6: The Schur-based viability kernel (safe) of Example 3.6.3 in transformed 2D subspaces.  3.6.4  4D Unstable System (An Example for Section 3.4)  Consider an unstable system [51, Ex. 2.2.1] with   0  1  0  0       0.1023 0 −0.0085 0  , A= 0 0 0 1   −0.0153 ε 0.0993 0   −3  B = 10  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  −1  −0.3426  0.0354 −0.6988    −0.0000 −0.2912 AT =   0 0  0  0  0.1399      0         0.9481 −0.0000  0  , T −1 B=10−3 ×   −0.7621 . 0.3150 −0.0135     0.0003 0.3188 0.3426  A target (unsafe) set K is chosen such that Zτ = {z ∈ R4 |  √  zTz ≤  0.2, z = T −1 x, 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 subsystem, to obtain an over-approximation of the full-order system, we only compute the over-approximation of the reachable tube in its stable subspace. 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 subsystems 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 resulting 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 decomposition 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 timedependent formulation of the disturbance to the upper subsystem and performing 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 Problem 2.1 draws upon the so-called Riccati transformation—a two-stage coordinate 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 standard Riccati transformation) which, in addition to parameterizing the unidirectional coupling between the subsystems, makes one of the subsystems ETUC. In the new coordinate space reachability analysis can then be performed in lower dimensions for each subsystem separately. The intersection of back projections of the lower dimensional reachable tubes is an overapproximation of the actual reachable tube in the transformed coordinate 1  A 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 A11 A12 B1  S=  (4.1)  A21 A22 B2  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 nonsingular transformation matrices T1 = T2 =  Ik  0  −L In−k Ik  M  0  In−k  ∈ Rn×n ,  (4.2)  ∈ Rn×n .  (4.3)  With L ∈ R(n−k)×k and M ∈ Rk×(n−k) that satisfy (NARE:) (Sylvester:)  R(L) := LA11 − A22 L − LA12 L + A21 = 0,  (4.4)  S (M ) := A11 − A12 L M − M A22 + LA12 + A12 = 0,  (4.5)  the transformed system is  A11 − A12 L A12 B1 , S = T1−1 (S) =  ✯0 ✟ ✟ R(L) A22 + LA12 LB1 + B2 ✟   S = T2−1 (S ) =  A11 − A12 L 0  ✿0 ✘✘ S✘ (M ) (I − M L)B1 − M B2 ✘  A22 + LA12  LB1 + B2  (4.6)  . (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 S1 =  A11 − A12 L (I − M L)B1 − M B2  S2 =  A22 + LA12 LB1 + B2  ,  .  (4.8) (4.9)  Suppose the control input is disjoint across these dynamically decoupled candidate 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. Parallelization of reachability calculations in each subspace could further reduce the computational time. Assumption 4.1. U =  2 i=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 evolution is tightly paired through common input. Difficulty arises, for example, when in the reachability computation a control value deemed op53  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 system. 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 (B2T ) ⊆ C (B1T ) 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 := −B2 B1† + Z − ZB1 B1† , 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 − B1 B1† A12 Z(B1 B1† − I)  (4.11)  + A22 − B2 B1† A12 Z(B1 B1† − I), where Π = −(B1 B1† − I) A11 + A12 B2 B1† ,  (4.12)  Γ = A22 B2 B1† + A21 − B2 B1† A12 B2 B1† + A11 .  (4.13)  To eliminate the non-invertible term (B1 B1† −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 − B1 B1† A12 Z + A22 − B2 B1† A12 Z  (4.14)  and δ ∈ R\{−1, 0} a finite (but possibly large) parameter such that B1 B1† −  (δ + 1)I is nonsingular:  R(Z) = ZΠ + Γ + Z A12 − B1 B1† A12 Z(B1 B1† − I) + A22 − B2 B1† A12 Z(B1 B1† − I)  (4.15)  . = ZΠ + Γ + F (Z)(B1 B1† − I) = δF (Z).  (4.16)  Simple algebraic manipulation and post-multiplication of R(Z)−δF (Z) = 0 by B1 B1† − (δ + 1)I  −1  then results in a NARE in the variable Z:  R1 (Z) := Z A˜11 − A˜22 Z − Z A˜12 Z + A˜21 = 0 −1 ˜ with A˜11 = Π B1 B1† − (δ + 1)I , A21 = Γ B1 B1† − (δ + 1)I B1 B1† A12 − A12 , and A˜22 = B2 B1† A12 − A22 .  (4.17) −1  , A˜12 =  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 − A22 L − LA12 L + 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 {−B2 B1† } 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 = T1−1 (S) =  A11 − A12 L  A12  B1  LA11 − A22 L − LA12 L + A21 A22 + LA12 LB1 + B2  (4.19)  =  A11 − A12 L δF (Z)  A12  B1  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 research [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 B2 B1† A12 − A22 is invertible. Define initial values as Z0 := B2 B1† A12 − A22 A0 := Π  B1 B1†  −1  − (δ + 1)I  Γ B1 B1† − (δ + 1)I −1  −  B1 B1† A12  −1  ,  (4.21)  − A12 Z0 .  (4.22)  To find Z we look for D := Z − Z0  (4.23)  by solving R1 (D) := DA0 − B2 B1† A12 −A22 + Z0 B1 B1† A12 − A12  D  − D B1 B1† A12 − A12 D + Z0 A0 = 0.  (4.24)  Lemma 4.2 ([62, Lem. 1]). Suppose B2 B1† A12 − A22 is nonsingular. If B2 B1† A12 − A22  −1  ≤  1 3  A0 +  B1 B1† A12  (4.25) − A12  Z0  then (4.24) has a unique real root D that satisfies 0≤ D ≤  2 A0 A0 +  Z0  B1 B1† A12  − A12  (4.26)  Z0  and is the fixed-point solution of the contraction Dk+1 = P1 (Dk ) given by P1 (Dk ) := B2 B1† A12 − A22  −1  Z0 A0 + Dk A0  − Z0 B1 B1† A12 − A12 Dk − Dk B1 B1† A12 − 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  B2 B1† A12 − A22  −1  A0 + B1 B1† A12 − A12  k  Z0  (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 − A12 L 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 = T2−1 (S ) =  A11 − A12 L − M δF (Z) δF (Z)  ✿0 ✘✘ R✘ 2 (M ) ✘  A22 + LA12 + δF (Z)M  B1 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 − A12 L be invertible and the initial  values be defined as  M0 := − A11 − A12 L  −1  A12 ,  (4.31)  N0 := A22 + LA12 + δF (Z)M0 .  (4.32)  We seek M by forming J := M − M0  (4.33)  and solving R2 (J) := JN0 − A11 −A12 L−δM0 F (Z) J +δJF (Z)J +M0 N0 = 0. (4.34) Lemma 4.3 ([62, Lem. 1]). Suppose A11 − A12 L is nonsingular. If A11 − A12 L  −1  ≤  1 3  (4.35)  N0 + δF (Z) M0  then (4.34) has a unique real root J that satisfies 0≤ J ≤  N0  2 N0 M0 + δF (Z) M0  (4.36)  and is the fixed-point solution of the contraction Jk+1 = P2 (Jk ) given by P2 (Jk ) := A11 − A12 L  −1  M0 N0 + Jk N0 + δM0 F (Z)Jk + δJk F (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 − A12 L  −1  k  N0 + δF (Z)  and decreases as δF (Z) , A22 , and  M0  A11 − A12 L  −1  (4.38) 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 fixedpoint 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 b |δ + 1|  ∀δ ∈ R\{−1, 0},  (4.39)  where the constants a and b are independent of δ and are determined by a := α(b/β)2 , b := 3 B1 B1† γβ, γ := A12 −  B1 B1† A12  , and β := A22 −  Γ  B2 B1† A12  .  A22 − B2 B1† A12  −1  , α :=  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 defeats the purpose of δ-parameterization of the unidirectional coupling term, since it can be shown that as δ grows, δF (Z(δ)) approaches a problemdependent 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 associated with the δ-parameterization in the modified Riccati approach and  instead enforces a trivial solution L = −B2 B1† . 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 f (δ) := δF (Z(δ))  δ∈R\{−1,0}  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 computed. 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 behavior similar to that of an absolute value proper rational function (over a  61  4.3. Recursive Decomposition  ||δF (Z(δ))|| ||Γ|| approximation  infeasible region  f(δ), fˆ(δ)  4  3.5  3  2.5  2 −1e6  −1e4  −1e2  1e2  1e4  1e6  δ  Figure 4.1: The worst-case unidirectional coupling f (δ) = δF (Z(δ)) (×’s) + 0.55| + 1.82 (dashed) and its approximation fˆ(δ) = | −27.65 δ 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 c0 fˆ(δ) = 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, consider 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 conservative 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 decomposition 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 −1 x, T = T1 T2 , the subsystem dynamics are governed by z˙1 = A11 − A12 L − δM F (Z) z1 + B1 u,  (4.41)  z˙2 = A22 + LA12 + δF (Z)M z2 + B2 u + δ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 −1 K 2:  Zτi ← ProjSi (Zτ ), ∀i ∈ {1, 2}  project onto ith subspace  For upper subsystem: 3:  1 1 Z[0,τ ] ← Reach[0,τ ] (Zτ , U)  For lower subsystem: 4:  Treat δF (Z)z1 as disturbance (existentially quantified)  5:  ζ ← z1 ≡ supv∈Z 1  [0,τ ]  6: 7: 8:  v  Compute upper-bound δF (Z)z1 ≤ δF (Z) ζ consrv.  2 2 Z[0,τ ] ←− Reach[0,τ ] (Zτ , B( δF (Z) ζ))  1 , Z2 ) return(Z[0,τ ] [0,τ ]  When the standard Riccati transformation is used to obtain the sub2 systems, steps 4–7 of Algorithm 4.1 are simply replaced with Z[0,τ ] ←  Reach[0,τ ] (Zτ2 , U).  Notice that if the error due to projections can be ignored, reachability in  1 −1 Reach the upper subspace is “exact” in the sense that Z[0,τ ] ≡ ProjS1 (T [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 subspace. 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,τ ] (Zτ2 , B( δF (Z) ζ)) ⊇ ProjS2 (T −1 Reach[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 directly proportional to supv∈Z 1  [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 N −1 2 Z[0,τ ]  = i=0  2 Z[iq, (i+1)q]  (4.44)  64  4.4. Reachability in Lower Dimensions and, by the semi-group property, 2 2 Z[iq, (i+1)q] = Reach[0,q] Z[(i+1)q, (i+2)q] , B( δF (Z) ζN −1−i ) ,  (4.45)  . 2 where Z[τ,+∞) = Zτ2 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 −1 sequence {ζi }N i=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 fullorder system in X can be obtained using Lemma 3.1 as in the Schur-based  case, i.e. via  1 2 T (Z[0,τ ] × S2 ) ∩ (S1 × Z[0,τ ] ) ⊇ Reach[0,τ ] (K, U).  4.4.1  (4.46)  Formulating an Upper-Bound on Conservatism of 2 Z[0,τ ]  2 2 Consider the computed reachable tube Z[0,τ ] := Reach[0,τ ] (Zτ , B( δF (Z) ζ))  in the lower subspace. Take z¯2 ∈ Zτ2 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 − 2 z2 ∈ Z[0,τ ]  τ  e−(τ −r)Ω ϑ(r)dr,  (4.47)  0  (4.48)  65  4.4. Reachability in Lower Dimensions with Ω := A22 + LA12 + δF (Z)M . Therefore, as in [36], z2 − e−τ Ω z¯2 ≤ =  τ  e(τ −r)  0 eτ Ω  Ω  Ω  −1  δF (Z) ζdr  (4.49)  δF (Z) ζ.  (4.50)  Notice that, eτ  Ω  Ω  −1  M  δF (Z) ζ = ≤  lim  M →∞  i=1 M  lim  M →∞  i=1  τi Ω i!  i−1  δF (Z) ζ  √ τ i σ(Ω) nk i!  (4.51)  i−1  δF (Z) ζ  =: η[0,τ ]  (4.52) (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  2  Z[0,τ ] ⊆  t∈[0,τ ]   e−Ωt Zτ2  ⊕ B(η[0,τ ] ).  (4.54)  2 The right-hand side of (4.54) provides an upper-bound on how much Z[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∈Z 1  [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 2  Z [0,τ ] :=  N −1     i=0  t∈[iq,(i+1)q]      2  ⊕ B(η[iq,(i+1)q] ) . (4.55) e−Ωt Z[(i+1)q,(i+2)q]  As in [36], it can be shown that the quality of this upper-bound is good in 66  4.4. Reachability in Lower Dimensions 2  2 ,Z the Hausdorff distance, in that limq→0 distH (Z[0,τ [0,τ ] ) → 0. Determining ]  2 −1 Reach whether Z[0,τ ] itself is a close over-approximation to ProjS2 (T [0,τ ] (K, U))  is much more involved and remains an open problem. We expect that performing reachability over sub-intervals, choosing nk appropriately, and minimizing 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 modified 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 dimensional 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 over67  Magnitude of unidirectional coupling  4.4. Reachability in Lower Dimensions  4 n=4 n=6 n=8 n=10 n=16  3.5 3 2.5 2 1.5 1 0.5 0 0.8  0.9  1  1.1 1.2 1.3 Decomposition time (s)  1.4  1.5  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 dimensional 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 −1 I = ∅ ⇐⇒ R[0,τ ] ∩ I = ∅ =⇒ R[0,τ ] ∩ I = ∅. 68  4.4. Reachability in Lower Dimensions Proof. ∀x ∈ R[0,τ ] ∀y ∈ I (R[0,τ ] ∩ I) = ∅ ⇐⇒ x = y ⇐⇒ T −1 x =  T −1 y ⇐⇒ T −1 R[0,τ ] ∩T −1 I = ∅. By Lemma 3.1 we have that R[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: 2  Proposition 4.5. i=1  i −1 I) = ∅ =⇒ R Z[0,τ [0,τ ] ∩ I = ∅. ] ∩ ProjSi (T  1 −1 I) = ∅)∧(Z 2 −1 I) = ∅) =⇒ (Z 1 Proof. (Z[0,τ ] ∩ProjS1 (T [0,τ ] ∩ProjS2 (T [0,τ ] × 2 ) ∩ T −1 I = ∅ ⇐⇒ Z −1 I = ∅ =⇒ R S2 ) ∩ (S1 × Z[0,τ [0,τ ] ∩ T [0,τ ] ∩ I = ∅. ]  That is, if the computed Riccati-based reachable tubes for each subsystem in the new coordinates do not intersect the projections of the transformed 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 c  tube R[0,τ ] , let R[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 conc  struct a feedback controller as in [81] on the boundaries of R[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 x˙ = Ax + Bu with   1.5072  3.3984    5.0644 −2.6683 A=  0.1156 −0.1863  −0.0808 0.0229   0.1300 −0.0884  0.0227 0.1689 , 0.5686 0.2648  0.4915 0.5949    −0.7433   −2.2528  B= −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  K} where ·  2  2  ≤ 0.2, z = T −1 x, x ∈  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   1.5362  3.4622  0  0    4.9377 −2.7030 A =  −1.8075 −0.0193  1.2341 0.2918  0  0  0.5727 0.4858      , 0.2612   0.5964    −0.7433       −2.2528  . B =   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 Example 4.5.1. each dimension for τ = 3 s. The full-order reachable tube (10144.80 s computational time) is over-approximated by the Riccati-based reachable tube (3.98 s computational time, including calculation of δ ∗ , transformation matrices, the decomposition, and projections) as shown in Figure 4.3. The loss of accuracy (due to treating the unidirectional coupling as disturbance 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    0  1  0  0      0.3920 0 −0.0327 0 , A=  0 0 0 1   0.0560 0 0.2753 0    0      −0.0033 . B=   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 −1 x, x ∈ K}.  We seek to identify the set of states for which there exists a bounded control law that keeps the system trajectories contained in Zτc . The safety-  preserving control synthesized through LS provides a guarantee that the pendulums’ angular displacement will not exceed an infinity norm ball around their upright positions, despite control saturation. Similar “envelope protection” 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 unidirectional coupling determined by the solution L = −B2 B1† 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   0  0.9524  0  0    0.3920 0 A =  0 0.1429   0  0  0  0      , 0 1.0500   0.2800 0    0       −0.0033  . B =   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 fullorder set. 73  4.5. Numerical Examples  4.5.3  Arbitrary 6D System  Consider the system x˙ = Ax + Bu with   3.3155    0.6320  −0.1047  A=  0.0344  −0.0140  −0.0106 B=  0.7768  0.0028 −0.0094 −0.0097  2.4455  −1.4796 −2.3001 −0.0370 −0.3522  0.6578  0.0282  0.0091  0.1885  −0.0225 −0.0012  0.2746  0.0360  −0.0215 −0.0082  0.0814   0.0322 −0.0112  −0.0621 −0.0214  , 0.0518 0.3567  0.0866 0.1567  0.0887 0.1182  0.1469 −0.7988 −2.3854 −0.0054  0.2657  and u ∈ R2 ,    1.3260 −0.1623  2.4582 −0.3955 −0.1403 −0.1187  T  0.3601  ≤ 0.05. We decompose this system into two 3D sub-  u  systems using the modified Riccati transformation with δ ∗ ≈ −199 (see Appendix B.3). A non-convex target set K is chosen such that in the trans3 j=1 Pj ,  formed coordinates Zτi = Czi − bj ≤ 0} with   1  C T = 0 0  −1 0 0  0 1 0  0 −1 0  0 0 1   0  0, −1  ∀i ∈ {1, 2}, where Pj = {zi ∈ R3 |  bT 1 = 0.5 bT 2  = −0.1  bT 3 = 0.5  0.3 0.3 0.3  0.1 0.1 0.1  0.1 0.1 0.1  −0.2  0.4 ,  0.4  0.4 ,  0.4  −0.2 ,  and z = [z1T , z2T ]T = T −1 x, 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 fullorder 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 solutions 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 algebraic 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 = B = [ 0.9600 0.1372 ]. Aric =  −2.0228 −0.3695  0.9732 0.0893  ,  Applying the modified Riccati transformation results in  −1.8360 0 −0.0875 −0.0975  , Bric = [ 0.9600 0 ]. Notice that the “fast” eigenvalue  λ1 = −1.8360 is assigned to the controlled subsystem and that the mag75  4.5. Numerical Examples  Success rate [%]  100  93.38%  88.08%  84.11%  79.47% 68.87%  50  0  4  6  8 Dimension  10  16  Figure 4.6: Fraction of a randomized test for which a Riccati solution existed. Success rate is shown as a percentage of 150 randomly generated two-time-scale systems for each dimension n = 4, 6, 8, 10, 16 with nk = n2 . 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 respectively. 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 Asch =  −0.0975 −0.1276 0 −1.8360  , Bsch =  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 Riccatibased decomposition is the superior approach in this example. 76  4.5. Numerical Examples 1.5  Schur  1 Riccati  0.5 0 −0.5 −1 −1.5 −1.5  −1  −0.5  0  0.5  1  1.5  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 decomposition 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 existentially quantified. This is done by replacing Step 3 in Algorithm 4.1 with 1 1 2 2 Z[0,τ ] ← Reach[0,τ ] (Zτ , U) (and Steps 4–7 with Z[0,τ ] ← Reach[0,τ ] (Zτ , U) if  the standard Riccati transformation is used to obtain the subsystems). As expected, we have ProjS1 (T −1 Reach[0,τ ] (K, U)) ≡ Reach[0,τ ] (Zτ1 , U),  (4.57)  ProjS2 (T −1 Reach[0,τ ] (K, U)) ⊆ Reach[0,τ ] (Zτ2 , B( δF (Z) ζ)).  (4.58)  and  77  4.5. Numerical Examples The over-approximation of the actual maximal reachable tube of the fullorder system can be obtained similarly to Lemma 3.1 as 2 1 T (Z[0,τ ] × S2 ) ∩ (S1 × Z[0,τ ] ) ⊇ Reach[0,τ ] (K, U).  (4.59)  While Algorithm 4.1 with the above modifications ensures that the fullorder (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 eigenvalues to the lower subsystem. Such an eigenvalue allocation is advantageous when computing the minimal reachable tube. In maximal reachability analysis, 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 supremum 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 approximated 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 x˙ = Ax + Bu, u ∈ R, |u| ≤ 0.1, with   0.5990 −0.2333 −0.0244  0.0121      −0.4553 −0.2107 −0.0076 −0.0041 ,  A= 0.0128 −0.1749 −0.1768 −0.0091  0.0309 −0.0229 −0.0105 −0.0777     0.6846    2.4813  . B=   0.3583 −1.5195 78  2  2  1.5  1.5  1  1  0.5  0.5 z4  z2  4.6. Summary and Further Discussions  0 −0.5  0 −0.5  −1  −1  −1.5  −1.5  −2 −2  −1  0 z1  1  (a) (z1 , z2 ) subspace  2  −2 −2  −1  0 z3  1  2  (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 decomposition (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 compares 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 lowerdimensional 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 reachability 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 subsystems, 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 minimal 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 reachable 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 becoming 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 transformation matrix. Finally, it is worth emphasizing that the proposed modified Riccati transformation 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 complexity 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 1  A 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 possibly 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 complexity.  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 x(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 sequence 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 }. We under-approximate K by a distance M P  (5.2)  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 t  x(t) − x(ti ) ≤  ti  x(τ ˙ ) 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 ) ∩ Reachtk −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 x˙ = 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 process, we reach an underapproximation K0 (P ) of the viability kernel.  Figure 5.1: Iteratively constructing V iab[0,τ ] (K, U).  an  under-approximation  of  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 x(τ ) M P  d ˚ K x(·)  ∂K  x0  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 ◦  ◦  x˙ = 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 → 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 ) ∈ Reachtk −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 ) ∈ Reacht1 −t0 (K1 (P ), U). So ◦  x0 ∈ K0 (P ). Hence V iab[0,τ ] (K, U) ⊆ 2  P ∈P([0,τ ]) K0 (P ).  See 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 results presented in Section 5.1.1. Define a sequence of sets recursively as  Kn = K,  (5.11a)  Kk−1 = K ∩ Reach1 (Kk , U) for k ∈ {1, . . . , n},  (5.11b)  where τ = n and Reach1 (·, ·) is the unit time-step maximal reachable set. Proposition 5.3. Let K0 be the final set obtained from the recurrence relation (5.11). Then, V iab[0,τ ] (K, U) = K0 .  (5.12)  Proof. Without loss of generality we assume that the time variable t is integer 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 ∈ Reach1 (Kk , U) for all k ∈ {1, . . . , n}. In particular, for k = 1 we have x0 := x(0) ∈ Reach1 (K1 , U). Thus, x0 ∈ K ∩ Reach1 (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 discretetime 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 available Lagrangian methods yield an (under- and/or over-) approximation of the maximal reachable set. The viability kernel should not be overapproximated 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 computation 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 reachability techniques that are based on polytopes (e.g. [76]) do not, in general, scale well with the dimension of the state. Moreover, the scalable reachability 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 dynamics 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, respectively. An ellipsoid in Rn is defined as E(q, Q) := x ∈ Rn | (x − q), Q−1 (x − q) ≤ 1 with center q ∈ Rn and shape matrix Rn×n  Q = QT  (5.14) 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 computations 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 [ δ]  Reachδ  δ  := (δ) ∈ M. We have that  (Kk (P ), U) = E(xc (0), X − (0)) ⊆  δ ∈M  E(xc (0), X − (0))  (5.17)  ⊆ Reachδ (Kk (P ), U). 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 continuoustime 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: K0 (P ) ← ∅ 4: while M = ∅ do 5: l← τ ∈M 6: k ← |P | 7: while k = 0 do 8: if Kk (P ) = ∅ then 9: K0 (P ) ← ∅ 10: break 11: end if [l] 12: G ← Reachtk −tk−1 (Kk (P ), U) Compute the ellipsoidal under-approximation of the maximal 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: K0∗ (P ) ← K0∗ (P ) ∪ K0 (P ) 17: M ← M\{l} 18: end while 19: return (K0∗ (P ))  91  5.2. Computational Algorithms Algorithm 5.2 Piecewise ellipsoidal approximation of V iab[0,τ ] (K, U) (discrete-time) 1: Kn ← K 2: K0∗ ← ∅ 3: while M = ∅ do 4: l← τ ∈M 5: k←n 6: while k = 0 do 7: if Kk = ∅ then 8: K0 ← ∅ 9: break 10: end if [l] 11: G ← Reach1 (Kk , U) 12: Kk−1 ← (Kn ∩ G) 13: k ←k−1 14: end while 15: K0∗ ← K0∗ ∪ K0 16: M ← M\{l} 17: end while 18: return (K0∗ ) Proposition 5.4. For a given partition P ∈ P([0, τ ]), let K0∗ (P ) be the set generated by Algorithm 5.1. Then,  K0∗ (P ) ⊆ V iab[0,τ ] (K, U).  (5.18)  Proof. Let K0 (P ) denote the final set constructed recursively by (5.4). Also, [l]  for a fixed direction l, let K0 (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, [l]  K0 (P ) ⊆ K0 (P ). Therefore, [l] l∈M K0 (P )  [l] l∈M K0 (P )  ⊆ V iab[0,τ ] (K, U).  ⊆ K0 (P ). Thus, K0∗ (P ) =  Remark 5.3. A similar argument holds for the discrete-time case in Algorithm 5.2, i.e. K0∗ ⊆ 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 := [l]  Reachtk −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 := [l]  Reach1 (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 underapproximated 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 | 1 2  x  2  ≤  1} with H = Q . Assume Y ∩ G = ∅ 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  log det H −1  (5.19a)  subject to  λi > 0   1 − λi 0 (q − qi )T   λi I H  0   (5.19b)  H∈Rn×n , q∈Rn , λi ∈R  q − qi  H  0, i = 1, 2.  (5.19c)  Qi  Using the optimal values for H and q, we will have (Y ∩ G) = E(q, H T H).  (5.20)  Loss of Accuracy A set generated by Algorithms 5.1 or 5.2 could be an inaccurate approximation of V iab[0,τ ] (K, U), especially for large time horizons. The loss of accuracy is mainly attributed to the function  (·), the under-approximation of  the intersection at every iteration with its maximum volume inscribed ellip93  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 intermediate 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 x(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 :=  0 E(0, [ 0.25 0 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 = [ 00 10 ] 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-approximation K↓ (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  |P | = 13  0.5  0.5  0  |P | = 21  0  -0.5 -0.5 0.5  0 0.5 |P | = 55  -0.5 -0.5 0.5  0  0 0.5 |P | = 89  -0.5 -0.5 0.5  -0.5 -0.5 0.5  0 0.5 |P | = 233  -0.5 -0.5 0.5  0 0.5 |P | = 144  0  0  0  |P | = 34  0.5  0 0.5 |P | = 377  -0.5 -0.5  0  0.5  0  0 -0.5 -0.5  0  0.5  -0.5 -0.5  0  0.5  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 Algorithm 5.1 under the double integrator dynamics. A finer time interval partition results in better approximation.  Error (# of grid points)  450  8.7%  400 350 6.1%  300 250  4.2%  200  3.6% 3.0%  150 100  13  21  34  55  |P |  89  3.0%  144  2.6% 233  2.4% 377  Figure 5.4: Convergence plot of the error as a function of |P | for the doubleintegrator 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 piecewise 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 x(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 ellipsoids over a finite subset of directions M ⊂ {l ∈ Rn | l, l = 1} and choose the one with the smallest volume:  El+  E(σ, Σ) := minvol AK ⊕ BU  .  l∈M  (5.23)  El+  Here we have used ⊕ 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 3. Under-approximate the set K  E( P σ, P  2  Σ).  2 Σ)  (5.24)  by computing tight  enclosed ellipsoids over directions in M and choose the one with the  largest volume to obtain K↓ (P ):  K↓ (P ) = maxvol K l∈M  Here the operator  El−  El−  E( P σ, P  2  Σ) .  (5.25)  returns the tight enclosed ellipsoidal approxi96  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 x(t) ˙ = Ax(t)+Bδe (t),     −0.003 0.039 0 −0.322 0.010     −0.065 −0.319  −0.180 7.740 0 , B =   A=  0.020 −0.101 −0.429 −1.160 0      0 0 1 0 0 ˙ θ]T ∈ R4 comprised of deviations in aircraft velocity with state x = [u, v, θ, [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     1075.84 0 0 0       0   0 67.24 0 0       , K = E     0 42.7716 0   2.18  0 0 0 0 76.0384 0   0  represents an ellipsoidal (under-)approximation of the flight envelope. We require x(t) ∈ K ∀t ∈ [0, 2]. 3  crad = 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 transparent region. The green/dark piecewise ellipsoidal sets underapproximate 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 underapproximation 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  5  5  5  0  0  -20  0  x1  -10  20  0  -5  -5  -5 -10  x4  10  x3  10  x2  10  -20  0  x1  -10  20  5  5  5  0  0  0  x2  10  -10 -10  0  x1  20  0  -5  -5  -5 -10 -10  x4  10  x4  10  x3  10  -20  0  x2  10  -10 -10  0  x3  10  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 kernel (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 delivery 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, mathematical guarantees of safety of the system are likely to be required. The viability kernel can provide such guarantees. Let us consider the problem of computing the viability kernel for a constrained dynamical system that represents the pharmacological response of a patient under anesthesia subject to therapeutic 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), y(t) = Cx(t)  (5.26) (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) ∈  4 The 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:   0.9960  0  0  0    0.0080 0.9960 0 0  −0.0080 0.0080 0.9960 0  A=  0.0079 −0.0080 0.0080 0.9960  −0.0079 0.0079 −0.0080 0.0080  0.0079 −0.0079 0.0079 −0.0080 B= C=  0.0894 −0.0890 18.5000  8.2300  0  0    0  0  0  0  0  0       ,      0.9960  0  0.0080  0.9960  0.0886 −0.0883  0.0879 −0.0876  3.5300  3.7000  4.3400  T  ,  3.0700 .  The constraint set (desired clinical effect) is specified in terms of the pseudooccupancy level and the input is bounded above and below by hard physical constraints:   y(t) ∈ K := [0.1, 1], 0 u(t) ∈ U := [0, 0.8].  (5.28)  0  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 coordinate 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 ···        · · · 0 −1 ∗ y y − y z1    ..     .       0 w2    w2  z2   ..   ..   ..  . .. .. . . .   .  =  .  =:  ..  .           z  ..  w w . 1 6 6      6 0 ∗ ∗ y y z7 ··· 0 1  (5.29)  Let u(t) = uss be the steady state control input needed for tracking a constant 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 reformulation, we deduct uss from the control set of the transformed system. The new constraints for the transformed, extended system z(t + 1) = Az(t) + Bu(t), y(t) = Cz(t), with A ∈ R7×7 , B ∈ R7×1 , C ∈ R1×7 , are as  follows:   z(t) ∈ K := (K − y ∗ ) × R6 , 0 u(t) ∈ U := U − u . 0  (5.31)  ss  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-approximate 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 approximate 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  z6  0  -20  -20 -20  0 20 z3  -20  -20  -20  0 20 z3  -20  0 20 z4  1  0  z6 0 20 z2  -20  0 20 z3  -20  0 20 z5  20  20  0  0 z1  -20 -20  0 20 z2  20  0  0  -20 -20  0 20 z2  20  20  0  -1  1  20  20  -20 -20  0 20 z2  0 z1  z4  0  -1  1  z6  -20  0 z1  20  -20  -20  -1  1  z5  z4  0  z5  z4 0 z1  20  20 z3  -1  1  0  -20  -20  z6  0 z1  0  z6  -1  0  -20  -20  -20  z5  0  z5  z3  z2  0  20  20  20  20  20  0  0  -20 -20  0 20 z4  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 piecewise 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 complexity 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 direction (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 control 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 differential 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 direction (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. x(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 ) ∩ Reachtk −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 Proposition 5.1: Since f is bounded on K, there exists a norm · and a real number 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 (nonanticipative) 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 concatenation 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 x˙ = 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 t  x(t) − x(tk ) ≤  tk  x(s) ˙ ds ≤ M (t − tk )  (6.4)  < M (tk+1 − tk ) ≤ M P . 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 Proposition 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 → 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 ) ∈ Reachtk −tk−1 (Kk (P ), U, V) for all k = 1, . . . , |P |, consider the tokenizations {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 ) ∈ Reacht1 −t0 (K1 (P ), U, V). So ◦  x0 ∈ K0 (P ). Hence Disc[0,τ ] (K, U, V) ⊆  P ∈P([0,τ ]) K0 (P ).  By approximating Disc[0,τ ] (K, U, V) using the sub-interval maximal reach-  able sets via recursion (6.2) we enable the use of scalable Lagrangian techniques for the computation of the discriminating kernel for high-dimensional systems. As we shall see, this will also allow us to compute the safetypreserving 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 recur109  6.2. Computational Algorithm & Control Synthesis sion (6.2). Corollary 6.2. For a fixed time partition P ∈ P([0, τ ]) suppose K0 (P ) = ∅.  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 Reachtk −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 x(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 (Algorithm 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 summarize the aforementioned algorithm used to under-approximate the discriminating 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 ∗[  ]  Kk−1τ =  [  for ∗[  ∗[  ]  τ (Kk (K|P | (P ) ∩ Reachtk −t k−1  τ]  (P ), U, V))  k ∈ {1, . . . , |P |}  (6.10)  ]  with K|P |τ (P ) = K|P | (P ) defined as in (6.2a) generates an ellipsoidal set ∗[  K0  τ]  (P ) such that ∗[  K0 τ ∈M  τ]  (P ) := K0∗ (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 K0∗ (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 K0∗ (P ) ⊆ K|P | (P ) ∩  [  ]  ∗[  Reacht1 τ (K1 τ ∈M  τ]  (P ), U, V)  (6.12)  ⊆ Disc[0,τ ] (K, U, V). Finally, note that similar to (6.6) the intermediate discriminating kernels are under-approximated via ∗[  τ ∈M  ]  ∗ Kk−1τ (P ) := Kk−1 (P ) ⊆ Disc[0,τ −tk−1 ] (K, U, V).  (6.13)  Remark 6.1. The under-approximation K|P | (P ) = K↓ (P ) may be constructed 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 El+  {AK, BU, GV} .  minvol l∈M  6.2.2  (6.14)  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 concatenating the sub-interval robust maximal reachability control laws according to Corollary 6.2. [  ]  ∗[  τ Suppose that in computing the reachable set Reachtk −t (Kk k−1  we can also compute [  ]  ∗[  τ Reachtk −t (Kk  τ]  (P ), U, V)  ∀t ∈ (tk−1 , tk ).  τ]  (P ), U, V)  (6.16)  That is, we compute the entire reachable tube over the kth sub-interval. For −[ τ ] ,k (t − tk−1 ) de[ τ] ∗[ ] note the center and the shape matrix of the ellipsoid Reachtk −t (Kk τ (P ), U, V) ∗[ ] c[ ] −[ ] at time t ∈ [tk−1 , tk ] with Kk τ (P ) = E(xk τ (tk − tk−1 ), X ,k τ (tk − tk−1 ))  fixed  c[  τ  ∈ M and k ∈ {1, . . . , |P |} let xk  and the terminal direction R(t, k) :=  τ ∈M  τ]  τ.  Define a shorthand notation  R[  τ]  (t, k) [  ]  ∗[  τ Reachtk −t (Kk  := τ ∈M  c[  = τ ∈M  (t − tk−1 ) and X  E(xk  τ]  τ]  (P ), U, V)  (t − tk−1 ), X  −[ τ ] ,k (t  (6.17)  − tk−1 )),  and suppose R(0, 1) = ∅. 112  6.2. Computational Algorithm & Control Synthesis  ∃ τ ∈ M s.t. x0 ∈ R[ τ ](0, 1) ∃ ¯τ = γ s.t. ◦ ¯ x(t) ∈ R[ τ ](t, k)  γ :=  τ  qperf  ◦  x(t) ∈ R(t, k)  uˆ(x(t), t) = ψBU (l[γ])  uˆ(x(t), t) ∈ U ◦  x(t) ∈ R[γ](t, k) γ := ¯τ  qsafe  ◦  ◦  x(t) ∈ R[γ](t, k)  x(t) ∈ R[γ](t, k)  γ := ¯τ  ∃ ¯τ = γ s.t. ◦ ¯ x(t) ∈ R[ τ ](t, k)  Figure 6.1: The graph of the hybrid automaton H representing the safetypreserving 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 (“performance” mode) and with qsafe representing the case in which the controller 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 direction 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 different 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 , ¯τ = γ}. The domains and the guards are chosen in terms of the interior 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 ∃ ¯τ = γ 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 by   U (q , γ, x(t), t) = U; fb perf U (q , γ, x(t), t) = ψ (l[γ] ), BU fb safe  (6.19)  where ψBU : Rn → U, l[γ] → µ − U B T l[γ] l[γ] , BU B T l[γ]  −1/2  (6.20)  is chosen so that BψBU (l[γ] ) is the support vector of the set BU = E(Bµ, BU B T ) ⊆ X in the direction −l[γ] ∈ Rn with l[γ] = l[γ] (x(t), t) = (X  −[γ] ,k (t  c[γ]  − tk−1 ))−1 (x(t) − xk (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 ∈ K0∗ (P ) where K0∗ (P ) is the piecewise ellipsoidal set obtained through (6.10)–(6.11), the feedback policy u(t) = u ˆ(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.  ∗ (P ) (which, as we shall see, First, note that if for every k, x(tk−1 ) ∈ Kk−1  ∗ (P ) ⊆ R(t is indeed the case) then x(tk−1 ) ∈ R(tk−1 , k) since Kk−1 k−1 , k).  Fix γ =  τ  ∈ M and k and let x(tk−1 ) ∈ R[γ] (tk−1 , k). Define a continuously [γ]  differentiable value function Vk : X × Tk → R such that [γ]  Vk (x(t), t) = dist2 (x(t), R[γ] (t, k)). [γ]  ◦  (6.23) [γ]  Notice that Vk (x(t), t) = 0 for x(t) ∈ R[γ] (t, k) and Vk (x(t), t) ≥ 0 for ◦  x(t) ∈ R[γ] (t, k). We use the convention  [γ]  R[γ] (tk−1 , k) ≡ x ∈ X | Vk (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) ∈ 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 Vk (x, t) along the system (6.9) yields d [γ] V (x(t), t) = 2 dist(x(t), R[γ] (t, k)) dt k d × dist(x(t), R[γ] (t, k)) dt  (6.25)  with ∂ d dist(x(t), R[γ] (t, k)) = dist(x(t), R[γ] (t, k)) dt ∂t ∂ + dist(x(t), R[γ] (t, k)), x(t) ˙ ∂x = l[γ] , Bu(t) + Gv(t) + ρBU (−l[γ] ) − l[γ] , Gv(t)  (6.26)  (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 dist(x(t), R[γ] (t, k)) = 0 dt d dist(x(t), R[γ] (t, k)) ≥ 0 dt  for u(t) = ψBU (l[γ] );  (6.29)  for any u(t) ∈ U.  (6.30)  When the disturbance does not play a worst-case game, (6.27) is replaced by  d dt  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  d dt  ∀v(t) = ρ[u](t) ∈ V.  dist(x(t), R[γ] (t, k)) ≤ 0 for u(t) = ψBU (l[γ] ) and ◦  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) = u ˆ(x(t), t) ∈ Ufb (q, γ, x(t), t)  we have  d [γ] V (x, t) dt k  (6.9)  ≤0  ∀(q, γ, x, t, v) ∈ Q × M × X × Tk × V.  (6.31)  This yields t tk−1  d [γ] V (x(s), s)ds ds k [γ]  [γ]  = Vk (x(t), t) − Vk (x(tk−1 ), tk−1 ) ≤ 0  (6.32)  which in turn implies (via (6.24)) [γ]  [γ]  Vk (x(t), t) ≤ Vk (x(tk−1 ), tk−1 ) ≤ 0  ∀t ∈ Tk .  (6.33)  Therefore, for every solution x of the differential inclusion x(t) ˙ ∈ Ax(t) ⊕ BUfb (q, γ, x(t), t) ⊕ GV we will have [γ]  ∗[γ]  x(t) ∈ R[γ] (t, k) = Reachtk −t (Kk (P ), U, V) ∀t ∈ Tk , ∗[γ]  x(tk ) ∈ Kk (P ).  (6.34) (6.35)  Notice that according to (6.13) ∗[γ]  Kk (P ) ⊆  ∗[γ]  Kk (P ) := Kk∗ (P ) γ∈M  (6.36)  ⊆ Disc[0,τ −tk ] (K, U, V) ⊆ K. 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 |} [γ]  ∗[γ]  Reachtk −t (Kk (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 ∈ u,ρ[u]  R(0, 1) we have xx0  6.2.3  (t) ∈ K ∀t ∈ [0, τ ] for every v(t) = ρ[u](t) ∈ V.  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 ) ∈ Kk∗ (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 Reachtk −t (Kk  τ]  (P ), U, V) at every time instant t, any one of these  ellipsoids can be used for the computation of l[ corresponding optimal control law ψBU  (l[  τ]  τ]  (x(t), t) in (6.21) and its  (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 Kk  τ]  (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 K0∗ (P ): Synthesis for x0 ∈ K0∗ (P ) When the initial condition x0 lies inside of the true discriminating kernel but outside of its piecewise ellipsoidal under-approximation, i.e. x0 ∈ Disc[0,τ ] (K, U, V)\K0∗ (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[  τ]  modified to be the direction that corresponds to any [  τ]  ∗[ dist(x(t), Reachtk −t (Kk  τ]  (x(t), t) in (6.21) is τ  ∈ M for which  (P ), U, V)) = min. We will demonstrate this us-  ing an example in Section 6.3.  Finally, note that for x0 ∈ 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 ∈ 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 controller returns Ufb (qsafe , γ, x(t), t) = ψBU (l[γ] ) and only this specific safetypreserving 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 between 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 frequency 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 ◦  ◦  −[γ] ,k (t − tk−1 ))  in qperf we  u(t) = (1 − βα [φ[γ] ](x(t), t))uperf (t) + βα [φ[γ] ](x(t), t)usafe (t),  (6.38)  c[γ]  and a given domain R[γ] (t, k) = E(xk (t − tk−1 ), X shall choose an input such that  with  βα [φ[γ] ](x(t), t) :=     1    1  1−α    0  (φ[γ] (x(t), t) − α)  if φ[γ] (x(t), t) ≥ 1;  if α ≤ φ[γ] (x(t), t) ≤ 1; if φ[γ] (x(t), t) ≤ α,  (6.39)  where α ∈ [0, 1) is a design parameter (Figure 6.2) and φ[γ] : Σe → R+ , (x(t), t) →  c[γ] (x(t) − xk (t − tk−1 )), −[γ] (X ,k (t − tk−1 ))−1 (x(t)  (6.40) (6.41) −  c[γ] xk (t  − tk−1 )) . ◦  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  0  α  1  φ  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 input 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 kernel via recursion (6.10) only the final intermediate maximal reachable sets [  ]  ∗[  τ Reachtk −t (Kk k−1  τ]  (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 [  ]  ∗[  τ Reachtk −t (Kk  τ]  t∈Tk  (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 (frequent switches between the two modes) or Zeno behavior (infinite switches in finite time). While such undesirable phenomena may be avoided by imposing 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 numerical 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 addition, 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 x˙ =  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. Figure 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 x˙ =  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  x0,1 0.2  x2  x0,4 0  x0,2  −0.2  x0,3 −0.4 −0.6 −0.6  −0.4  −0.2  0  x1  0.2  0.4  0.6  Figure 6.3: Closed-loop trajectories in the phase-plane for four sample initial 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 approximation 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 finite 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 trajectories can leave the kernel since it is finite-horizon and hence not invariant; however, they do not violate the constraints. For x0,4 ∈ 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 trajectory eventually leaves K.  124  6.3. Numerical Examples  u  0.1 0 −0.1  (x0,1 ) 0  0.5  1  1.5  2  u  0.1 0 −0.1  (x0,2 ) 0  0.5  1  1.5  2  0.1 0 −0.1  u  (x0,3 ) 0  0.5  1  1.5  2  0.1 0 −0.1  u  (x0,4 ) 0  0.5  1  1.5  2  t Figure 6.4: The corresponding safety-preserving feedback policy for Example 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  x0,1 0.2  x2  x0,4 0  x0,2  −0.2  x0,3 −0.4 −0.6 −0.6  −0.4  −0.2  0  x1  0.2  0.4  0.6  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  u  0.1 0 −0.1  (x0,1 ) 0  0.5  1  1.5  2  u  0.1 0 −0.1  (x0,2 ) 0  0.5  1  1.5  2  0.1 0 −0.1  u  (x0,3 ) 0  0.5  1  1.5  2  0.1 0 −0.1  u  (x0,4 ) 0  0.5  1  1.5  2  t 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  x0,2 0.2  x2  x0,3 0  −0.2  x0,1  −0.4 −0.6 −0.6  −0.4  −0.2  0  x1  0.2  0.4  0.6  Figure 6.7: Closed-loop trajectories in the phase-plane for three sample initial 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 disturbance.  128  6.3. Numerical Examples  u  0.2  (x0,1 )  0  −0.2  0  0.5  1  1.5  2  u  0.2  (x0,2 )  0  −0.2  0  0.5  1  1.5  2  u  0.2  (x0,3 )  0  −0.2  0  0.5  1  1.5  2  t Figure 6.8: The corresponding safety-preserving feedback policy for Example 6.3.2 for each initial condition x0,i .  v  0.1  (x0,1 )  0  −0.1 0  0.5  1  1.5  2  v  0.1  (x0,2 )  0  −0.1 0  0.5  1  1.5  2  v  0.1  (x0,3 )  0  −0.1 0  0.5  1  1.5  2  t Figure 6.9: The simulated disturbance signal for Example 6.3.2 for each initial condition x0,i .  129  6.3. Numerical Examples k10 0.4436  k12 0.1140  k13 0.0419  k21 0.0550  k31 0.0033  V1 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 Section 5.3.2. Instead of the discrete-time Laguerre model of the patient, here we use a three-dimensional continuous-time pharmacokinetic/pharmacodynamic (PKPD) compartmental model whose state variables describe the concentration of propofol in each compartment of the body:     −(k10 + k12 + k13 ) k12 k13 1/V1     x˙ =  k21 −k21 0  x +  0  u. k31  0  −k31  0  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    3.5   x(t) ∈ K := E  5  , 5   6.25 0  25  0 0  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 T  the initial condition x0 = 4 5 8 and assume that no drug is being administered, 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 10  x3  x2  10  5  0  0  2  4  x1  5  0  6  0  2  x1  4  6  x3  10  5  0  0  5  x2  10  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 kernel (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 application 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  10  x3  x2  10  5  0  0  2  4  x1  5  0  6  0  2  x1  4  6  x3  10  5  0  0  5  10  x2  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 kernel (green/light) are also shown. Safety is preserved.  20  u  15 10 5 0 0  5  10  15 t  20  25  30  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  100 % uperf in u  80 60 40 20 0 0  5  10  15 t  20  25  30  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 component 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 correspond 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.  active mode  safe  perf 0  5  10  15 t  20  25  30  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  15  γ  10  5  0  5  10  15 t  20  25  30  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.  1  φ  0.8 0.6 0.4 0.2 0  0  5  10  15 t  20  25  30  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  10  x3  x2  10  5  0  0  2  4  x1  5  0  6  0  2  x1  4  6  x3  10  5  0  0  5  10  x2  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.  20  u  15 10 5 0 0  5  10  15 t  20  25  30  Figure 6.18: The safety-preserving feedback policy for Figure 6.17 for Example 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  active mode  safe  perf 0  5  10  15 t  20  25  30  Figure 6.19: The active mode of the hybrid automaton H at time t for Example 6.3.3 when the modified policy (6.38)–(6.41) is not employed.  15  γ  10  5  0  5  10  15 t  20  25  30  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  1  φ  0.8 0.6 0.4 0.2 0  0  5  10  15 t  20  25  30  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 x˙ = Ax + Bu with  −0.0226 −36.6170 −18.8970 −32.0900 3.2509   0.0001 −1.8997 0.9831 −0.0007 −0.1708   0.0123 11.7200 −2.6316 0.0009 −31.6040  A=  0 0 1.0000 0 0   0 0 0 0 −30.0000  0  B=  0  0  0 0 0 0 30  0  0 0 0 0  30  0  0  −0.7626     −0.0050   22.3960  ,  0   0  −30.0000  0  T  T  and state vector x = α˙ α θ˙ θ xδe xδc 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. T  The control input u = δe δc 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 ∞  (xT Qx + uT Ru)dt  J(u) =  (6.44)  0  with Q = I6 and R = 102 × 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.0552  0.5619  0.2367  0.3119 −0.2107  0.3588 −0.2107  0.1497  .  Since the control authority is constrained by the ellipsoid U =: E(ω, Ω), 1  These assumptions are purely for the algorithmic convenience of dealing with ellipsoids. 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  u lqr sat(ulqr ) := ω + Ωu  lqr  ulqr , Ωulqr  −1/2  if ulqr ∈ U;  (6.45)  if ulqr ∈ U.  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 disturbance/uncertainty. Consequently, the discriminating kernel is expressed in terms of robust maximal reachable sets. Owing to this new connection, scalable and efficient Lagrangian methods can now be used to correctly approximate the discriminating kernel. We also presented an extension 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  5 x3  x2  5 0  -5 -20  -10  x1  0  10  -5 -20  0  -10  x1  0  10  0  10  0 x2  5  0 x3  5  5 x4  x4  x1  0  -5 -5  5 0  -5 -5  -10  5 x3  x4  5  -5 -20  0  0 x2  5  0  -5 -5  Figure 6.22: 2D projections of the closed-loop trajectories with saturated LQR and no safety control for the first four states for Example 6.3.4. The initial condition x0 is marked by ‘×’. The constraint set K (red/dark) and a piecewise ellipsoidal underapproximation 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 x3  x2  5 0  -5 -5  0 x1  5  -5 -5  0  0 x1  5  0 x2  5  0 x3  5  5 x4  x4  5  0  -5 -5  5 0  -5 -5  0 x1  5 x3  x4  5  -5 -5  0  0 x2  5  0  -5 -5  Figure 6.23: 2D projections of the closed-loop trajectories with safetypreserving 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 piecewise 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 demonstrated 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. Reachability 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 approximate 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 higherdimensional systems. The first approach, based on structure decomposition, aims to facilitate the use of Eulerian methods on higher-dimensional, continuous-time, continuously-valued LTI systems (and by extension, hybrid 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 subspaces. It is shown that the reverse transformation of the intersection of the back-projection of these resulting sets correctly approximates the actual constructs. Within the framework of structure decomposition, we proposed two techniques: the Schur-based decomposition and the Riccati-based decomposition, each with its own merits. While the Schur-based decomposition is quite generic and applicable to most systems, the Riccati-based decomposition 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 methods, draws a connection between the viability kernel and the maximal reachable 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 underapproximated with polynomial complexity. Using the well-established ellipsoidal techniques for maximal reachability we then proposed a scalable algorithm that facilitates the computation of the viability kernel for highdimensional 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 disturbances or uncertainties. We demonstrated our techniques on several examples (up to 8D) including 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 decomposition 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 algorithm which was proposed based on this new connection using the ellipsoidal techniques for maximal reachability is capable of computing 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 scalable 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 particularly 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 natural (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 discriminating kernel case) may be used to formulate a safety-preserving controller based on the framework described in Chapter 6. Our proposed control algorithm 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 currently 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 converge to the infinite-horizon kernel. However, if the algorithm does converge for at least one terminal direction, then the current formulation 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 control 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 analysis. In Proc. IEEE International Symposium on Computer-Aided Control Systems Design, Munich, Germany, October 2006. [5] J.-P. Aubin. Viability Theory. Systems and Control: Foundations and Applications. Birkh¨ auser, Boston, MA, 1991. [6] J.-P. Aubin. Viability kernels and capture basins of sets under differential 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 autolander safety analysis through optimal control-based reach set computation. Journal of Guidance, Control, and Dynamics, 30(1):68–77, 2007. [9] C. B´en´e, L. Doyen, and D. Gabay. A viability analysis for a bioeconomic 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. Conference 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 Annals of the International Society of Dynamic Games, pages 177–247, Boston, MA, 1999. Birkh¨auser. [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 Approximate 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 Hybrid Systems: Computation and Control, LNCS 1386, pages 96–109. Springer, 1998. [23] A. N. Daryin, A. B. Kurzhanski, and I. V. Vostrikov. Reachability approaches 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 kernels 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 Discontinuous 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 Verification (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 diagonalization 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: Computation 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 approximate 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. Hespanha and A. Tiwari, editors, Hybrid Systems: Computation and Control, 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 Hopkins 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. Automated 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 systems using reduced-order models. In Proc. American Control Conference, 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. Tiwari, 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 Control 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 propofol 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ˇci´c, and S. Ramesh. Generating and analyzing symbolic traces of Simulink/Stateflow models. In A. Bouajjani 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 Hybrid Systems: Computation and Control, pages 55–63, Beijing, China, 2012. [56] 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; 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 Schurbased 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 continual 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. Fixedcomplexity piecewise ellipsoidal representation of the continual reachability 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´c. A Riccati equation for block-diagonalization of ill-conditioned systems. IEEE Transactions on Automatic Control, 20(6):812–814, 1975. [63] P. V. Kokotovi´c, 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: optimization and parallel computations. Optimization methods and software, 14(4):267–310, 2001. [65] B. H. Krogh and O. Stursberg. Efficient representation and computation 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¨ uhn. 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 techniques 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´alyi. Ellipsoidal Calculus for Estimation and Control. Birkh¨ auser, Boston, MA, 1996. [70] A. B. Kurzhanski and P. Varaiya. Ellipsoidal techniques for reachability analysis. In N. Lynch and B. Krogh, editors, Hybrid Systems: Computation and Control, LNCS 1790, pages 202–214, Berlin Heidelberg, 2000. Springer-Verlag. [71] A. B. Kurzhanski and P. Varaiya. Ellipsoidal techniques for reachability 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 reachability 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´c, and M. Morari. Multi-Parametric Toolbox (MPT). In R. Alur and G. J. Pappas, editors, Hybrid Systems: Computation and Control, LNCS 2993, pages 448–462, Berlin, Germany, 2004. Springer. [77] C. Le Guernic. Reachability analysis of hybrid systems with linear continuous dynamics. PhD thesis, Universit´e Grenoble 1 – Joseph Fourier, 2009. [78] C. Le Guernic and A. Girard. Reachability analysis of linear systems using support functions. Nonlinear Analysis: Hybrid Systems, 4(2):250–262, 2010. [79] J. Lygeros. On reachability and minimum cost optimal control. Automatica, 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. Constrained model predictive control: Stability and optimality. Automatica, 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 nonlinear systems with terminal integrators: a mixed implicit explicit formulation. 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. SpringerVerlag. [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 Automatic 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¸ca. Supervised 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. Kyriakopoulos. A viability approach for the stabilization of an underactuated 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´c. Set theoretic methods in model predictive control. Nonlinear Model Predictive Control, pages 41–54, 2009. [101] S. V. Rakovi´c 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 balancedtruncation model reduction. IEEE Transactions on Automatic Control, 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 Mathematics and Optimization, 29(2):187–209, Mar 1994. 157  Bibliography [107] K.-H. Shim and M. E. Sawan. Singularly perturbed unified time systems 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. Hofmockel, and B. Lampe. Automatic drug delivery in anesthesia: From the beginning until now. In Proc. Mediterranean Conf. Contr. Automation, Athens, Greece, 2007. [110] J. M. Siret, G. Michailesco, and P. Bertrand. Optimal approximation 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 transformation. SIAM Review, 29(1):91–113, 1987. [113] D. M. Stipanovi´c, I. Hwang, and C. J. Tomlin. Computation of an overapproximation 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˜ no, C. Ionescu, and R. De Keyser. NMPC for propofol drug dosing during anesthesia induction. In Nonlinear Model Predictive 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. Computational 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 modeling. 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 (B T ) ⊆ C (AT ).  Proof. By taking the transpose of both sides we have AT X T = B T . Let Y := X T . The equation AT Y = B T is known to have at least one solution w.r.t. Y if and only if the column space of B T is in the image (or range) of the linear operator AT , i.e. C (B T ) ⊆ 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  A.2  AT | B T  = rank AT .  (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) A (K) over the interval [0, t], t ∈ [0, τ ], for a fixed τ ∈ R+ . Denote by Reach[0,t]  160  A.2. Proof of Proposition 3.4 the backward reachable tube of its corresponding autonomous system x˙ = Ax.  (A.3)  The following inclusions hold: A K ⊆ Reach[0,t] (K, U) ⊆ Reach[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 =  p i=1 Ui ,  ui ∈ Ui = [U i , U i ], 0 ∈ Ui . Notice that  the trajectories of the autonomous system (A.3) approach those of the controlled 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 tube Reach[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)  H(x, ) = sup , Ax + Bu  (A.6) (A.7)  u∈U  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, ) =  T  Ax +  T  Bu∗ ,  u∗ = [u∗1 · · · u∗p ]T  (A.8)  161  A.2. Proof of Proposition 3.4 with u∗i  =     U   i  if  Tb i  < 0;  [U i , U i ] if    U if i  Tb i  = 0;  Tb i  >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. Therefore we have  0 ∇t φ(x, t) = | T Ax| −  if T Bu∗  T Ax  ≥−  T Bu∗ ;  T Bu∗  ≥ 0.  (A.10)  otherwise.  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, ) = T Ax.  Consequently, (A.10) reduces to  ∇t φ(x, t)  ξ=0   0 if T Ax ≥ 0; =: ∇t φA (x, t) = | T Ax| otherwise  (A.11)  where φA (·, ·) is to denote the implicit surface function whose zero sub-level  A (K) of (A.3). That is, set determines the backward reachable tube Reach[0,t] A (K) = {x | φA (x, τ − t) ≤ 0}. Reach[0,t]  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. T Ax  ≥ 0 as opposed to  value is increased (i.e. |  T Ax  T Ax|  ≥−  T Bu∗ ),  but also its maximum (positive)  as opposed to |  for all (x, t) ∈ Rn × [0, τ ] we have  T Ax|  −  T Bu∗ ).  ∇t φA (x, τ − t) ≥ ∇t φ(x, τ − t)  =⇒ φA (x, τ − t) ≤ φ(x, τ − t) ≤ φ(x, τ ) ≤ 0 A ⇐⇒ Reach[0,t] (K) ⊇ Reach[0,t] (K, U) ⊇ K.  Therefore,  (A.12) (A.13) (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. ξ = 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 A (K) ⊇ Reach[0,t] (K, U) ∀t Reach[0,t]  (A.15)  A (K) denotes the backward reachable tube of the autonomous where Reach[0,t]  system (A.3). Therefore, to prove Reach[0,t] (K, U) = K, ∀t ∈ [0, τ ], it is A (K) = K, ∀t ∈ [0, τ ]. sufficient to show that Reach[0,t]  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  = ∇x φ (x, t), SΛS  (A.16) −1  x  (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 Reach[0,t] (K)  = K.  (A.19) (A.20)  This concludes the proof.  163  A.3. Decomposed System Matrices for Examples 3.6.2 and 3.6.3 x2 A2  x1 A 1  x3 A3  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:   0.1527 −0.1511  0.0312  0       0.1853 −0.1536 −0.0247 0.0065   , Ad =   0 0 −0.3169 −7.5973   0 0 0.1028 −0.4331    0        0 .  Bd =   −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:   −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  Ad =   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   0 0     0 0     0 0     0   0 Bd =  .  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   Ad2    =    −1.3594 −1.2819  0  0       1.0769 −0.3660  0 0 ,  =  0 0 −0.9978 −2.8164   0 0 3.0041 −0.2943 −0.0052 0 0 0  0.0026  0.0215  0.5192     Bd1    0 0      0 0 ,  =  0 0   0 0 0  0    −0.0755 −0.1334 0 0      0  0.1262  0  , Bd =  . 2    −0.7479 −0.1877  1.3130 −0.0574   0.0339 −0.2137 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 +U CV )−1 = Y −1 − Y −1 U (C −1 + V Y −1 U )−1 V Y −1 , with Y = −(δ + 1)I, U = B1 , C = I, and V = B1† we have  −1  B1 B1† − (δ + 1)I  =−  1 1 I + B1 B1† . δ+1 δ  (B.1)  Using this, (4.14), (4.21), (4.23), (4.26), multiplicative and triangular inequalities, and the fact that B1 B1† ≥ 1 we obtain δF (Z(δ)) ≤ |δ| α( Z0 + D )2 + β( Z0 + D ) ≤ |δ| α  Z0 + +β  ≤ |δ| 9α Z0 ≤ |δ| 9αγ 2  2  Z0 +  2 A0 Z0 A0 + α Z 0  + 3β Z0  B1 B1† − (δ + 1)I + 3βγ  ≤ |δ| 9αγ 2  2  2 A0 Z 0 A0 + α Z0  1 δ+1  2  + 3βγ  −1 2  B1 B1† − (δ + 1)I 1+  1 δ  1 δ+1  2  B1 B1†  1+  1 δ  −1  2  B1 B1† 166  B.2. Condition Number of the Modified Riccati Transformation  =  1 |δ| + 1 |δ| |δ + 1|  2  a+  |δ| + 1 b |δ + 1|  ∀δ ∈ 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  δ→±∞  δ 1 Q1 (I + B1 B1† )P1 Q1 2 δ→±∞ (δ + 1) δ δ 1 1 × (I + B1 B1† ) + P2 Q1 (I + B1 B1† ) δ δ+1 δ  δF (Z(δ)) = lim  (B.2)  = 0 + P2 Q1 = Γ −1  B2 B1† A12 − A22  with Q1 :=  (B2 B1† A12 − A22 ).  B.2  (B.3)  Γ, P1 := (A12 − B1 B1† A12 ), and P2 :=  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 = T1 T2 is bounded by κ(T ) ≤ max 1 + µ, 1 + ν(1 + µ) · max 1 + ν, 1 + µ(1 + ν)  (B.4)  with constants µ :=  N0  2 N0 M0 + δF (Z)  M0  + M0 ,  ν := B2 B1† + (1 + B1 B1† )  2 A0 A0 +  Z0  B1 B1† A12  − A12  Z0  Proof. Let A = [Aij ] be any partitioned matrix. Then A ≤  + Z0 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 =  T −1 ≤  κ(T ) = T  1 L  I−M L −M L I  , we have 1+ M L L  M 1+ L M  From (4.10) we find L ≤ B2 B1† + Z  1 + B1 B1†  and (4.26) we know that Z is bounded above by Z0 , and similarly from (4.33) and (4.36), 2 N 0 M0 N0 + δF (Z) M0  + M0 . Therefore,  M  M 1  .  (B.5)  . But from (4.23)  2 A0 Z0 A0 + B1 B1† A12 −A12  Z0  +  is bounded above by  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:   3.3126  0.7676  2.4511    0.6223 −1.5072 −2.3105    −0.0989 −0.3285 0.6852 A =  0.0944 0.0042 −0.1299    −0.0540 −0.4170 −0.0461  −0.1474 0.1949 0.2443   0.1469  0  0  0        0 0 0  , 0.1880 0.0445 0.3528    0.2802 0.0888 0.1593   0.0848 0.0888 0.1197 0  0.2657  0  0       −0.7988 2.4582       −2.3854 −0.3955   . B =  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:   0.5912 −0.2477  0  0       −0.4583 −0.2017  0 0  A =  0.0197 0.1351 −0.1905 −0.1878  ,   −0.0825 −0.1479  0.0012 −0.0633    0.6846       2.4813   B =  0 .   0  169  Appendix C  Other Backward Reachability Constructs1 Some additional backward constructs formed under the constrained dynamical 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 trajectories emanating from those states reach K exactly at time t: Reacht (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 reachable 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 1  A 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 mentioned in Section 2.2: Proposition C.3 ([86]). Reach[0,τ ] (K, U) ⊇  t∈[0,τ ]  Reacht (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) = Proof. x0 ∈  t∈[0,τ ] Reacht (K, U)  K ⇐⇒ ∀u(·) ∈ U[0,τ ] ∀t ∈ [0, τ ]  t∈[0,τ ]  Reacht (K, U).  (C.7)  ⇐⇒ ∀t ∈ [0, τ ] ∀u(·) ∈ U[0,τ ] xux0 (t) ∈ xux0 (t) ∈ K ⇐⇒ x0 ∈ Inv[0,τ ] (K, U).  This can also be verified from (C.5) and (2.9) and the simple fact that Reacht (K, U) = (Reacht (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,τ ]  Reacht (K, U).  (C.8)  Proof. x0 ∈ Reachγ[0,τ ] (K, U) ⇐⇒ ∀t ∈ [0, τ ] ∃u(·) ∈ U[0,t] xux0 (t) ∈ K ⇐⇒ ∀t ∈ [0, τ ] x0 ∈ Reacht (K, U) ⇐⇒ x0 ∈  t∈[0,τ ] Reacht (K, U).  172  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0072951/manifest

Comment

Related Items