A Modified Riccati Transformation for Decentralized Computation of the Viability Kernel Under LTI DynamicsResearch supported by NSERC Discovery Grant #327387 (M. Oishi), NSERC Collaborative Health Research Project #CHRPJ-350866-08 (G. Dumont), and the Institute for Computing, Information and Cognitive Systems (ICICS). This work was mainly carried out at Electrical & Computer Engineering, University of British Columbia, Vancouver, BC V6T 1Z4, Canada.

A Modified Riccati Transformation for Decentralized Computation of the Viability Kernel Under LTI Dynamicsthanks: Research supported by NSERC Discovery Grant #327387 (M. Oishi), NSERC Collaborative Health Research Project #CHRPJ-350866-08 (G. Dumont), and the Institute for Computing, Information and Cognitive Systems (ICICS). This work was mainly carried out at Electrical & Computer Engineering, University of British Columbia, Vancouver, BC V6T 1Z4, Canada.

Shahab Kaynama  and Meeko Oishi S. Kaynama (kaynama@ece.ubc.ca, cor. author) is currently with Electrical Engineering & Computer Sciences, University of California at Berkeley, 337 Cory Hall, Berkeley, CA 94720, USA.M. Oishi (oishi@unm.edu) is with Electrical & Computer Engineering, University of New Mexico, MSC01 1100, 1 University of New Mexico, Albuquerque, NM 87131, USA.
(Preprint Submitted for Publication)

Computing the viability kernel is key in providing guarantees of safety and proving existence of safety-preserving controllers for constrained dynamical systems. Current numerical techniques that approximate this construct suffer from a complexity that is exponential in the dimension of the state. We study conditions under which a linear time-invariant (LTI) system can be suitably decomposed into lower-dimensional subsystems so as to admit a conservative computation of the viability kernel in a decentralized fashion in subspaces. We then present an isomorphism that imposes these desired conditions, particularly on two-time-scale systems. Decentralized computations are performed in the transformed coordinates, yielding a conservative approximation of the viability kernel in the original state space. Significant reduction of complexity can be achieved, allowing the previously inapplicable tools to be employed for treatment of higher-dimensional systems. We show the results on two examples including a 6D system.

1 Introduction

Constrained dynamical systems have received a tremendous amount of attention due to the presence of safety constraints and hard bounds that appear in many practical scenarios. Providing guarantees of constraint satisfaction and facilitating synthesis of constraint-satisfying controllers therefore is highly desirable, particularly in safety-critical applications. A class of safety-critical systems known as envelope protection problems is concerned with ensuring that the trajectories remain in a safe, bounded “envelope” (subset) of the state space for a given time horizon. Such problems arise in e.g. flight management systems [1, 2, 3, 4] where the safety constraints are defined as the aircraft’s aerodynamic envelope and consequently the system must ensure that certain combinations of states are avoided to prevent stalling or other undesirable behaviors. Other application domains include control of depth of anesthesia [5], aircraft autolanders [6], automated highway systems [7], control of under-actuated underwater vehicles [8], stockout prevention of storage systems in manufacturing processes [9], and management of a marine renewable resource [10], to name a few.

Viability theory [11, 12, 13] provides a set-valued perspective on the behavior of the trajectories inside a given set. Thus it is naturally suited to handle envelope protection problems. By duality, minimal reachability [14] is also capable of analyzing such problems by investigating the behavior of the trajectories outside of the envelope. For simplicity, in this paper we only focus on the constructs generated within the framework of viability theory. The viability kernel is the set of initial states for which there exists at least one trajectory of the input-constrained system that respects the state constraint for all time. It is shown in [12] and (by duality in [15]) that the viability kernel is the only construct that can be used to prove safety/viability of the system and to synthesize inputs that preserve this safety; cf. [16, Chap. 1–2] for more detail. In general an exact computation of the viability kernel is extremely difficult if not impossible. Instead, approximations of this set are computed. Such computations have historically been subject to Bellman’s “curse of dimensionality” [17]. The numerical algorithms that approximate the viability kernel and its associated control laws (e.g., [14, 18, 19, 20]), collectively referred to as Eulerian methods [15], rely on gridding the state space and therefore their computational complexity increases exponentially with the dimension of the state. This renders them impractical for systems of dimension higher than three or four.

This paper presents a part of our efforts to address the curse of dimensionality by enabling the use of Eulerian algorithms for higher-dimensional LTI systems (and by extension, hybrid systems with LTI dynamics). We decompose the structure of the system, applying Eulerian algorithms on each individual lower-dimensional subsystem in a decentralized fashion. Significant computational gains can be obtained, since instead of one costly centralized computation on the full-order system, multiple less expensive subsystem computations are performed. The results are then mapped back to the full-order space to obtain a conservative approximation (i.e. an under-approximation) of the viability kernel. The contribution of this paper is twofold: 1) We investigate various structures on system matrices that must be satisfied so that the behavior of the constrained system for envelope protection problems (with simply-connected, compact constraints) can be inferred conservatively from subspace decentralized analyses (Section 3). 2) We then present an isomorphism through which the desired structure is imposed on the system (albeit under certain conditions) to facilitate decentralized computations in the transformed space (Section 4). Numerical examples are provided in Section 5.

1.1 Related Work

Complexity reduction for viability and minimal reachability has been addressed by many researchers. A projection scheme in [21] based on Hamilton-Jacobi (HJ) partial differential equations (PDEs) over-approximates the projection of the true minimal reachable tube in lower dimensional subspaces, with the unmodeled dimensions treated as a disturbance. Similarly, [22] decomposes a full-order nonlinear system into either disjoint or overlapping subsystems and solves multiple HJ PDEs in lower dimensions. More recently, a mixed implicit-explicit HJ is presented in [23] for nonlinear systems whose state vector contains states that are integrators of other states. The complexity of this new formulation is linear in the number of integrator states, while still exponential in the dimension of the rest of the states. These techniques assume that the system itself presents a certain structure that can be exploited.

In [24], an approximate dynamic programming technique is presented that, although still grid-based, enables a more efficient computation of the viability kernel. The viability kernel (similarly to [25]) is expressed as the zero sublevel set of the value function of the corresponding optimal control problem. It is assumed that the value function, which is a viscosity solution of a HJB PDE, is differentiable everywhere on the constraint set. The PDE is then discretized and the resulting value function is numerically computed on a grid using a function approximator such as the -nearest neighbor algorithm. The error-bounded approximation is not conservative (it is an over-approximation) but converges to the true viability kernel in the limit as the number of grid points goes to infinity.

Another related approach is the search for a barrier certificate [26], a Lyapunov-like function that forms a separating hyper-surface between any two given sets and in the state space. If there exists a function non-positive on and positive on , and whose Lie derivative (along the vector field) is non-positive on its zero level set for all states and controls, then no trajectories will ever go from to . This technique can be adapted to analytically describe the boundary of the infinite-horizon viability kernel: A certificate must now be formulated such that at every state along its zero level set there exists a control that makes the Lie derivative non-positive. For systems with polynomial vector fields and semi-algebraic constraints, efficient techniques based on Sum of Squares can be used to find the barrier certificate.111This method cannot be used to formulate the finite-horizon viability kernel which may be useful when, for example, the infinite-horizon kernel is empty, or when safety is to be verified/enforced over a finite time interval. Moreover, there are no guarantees that a barrier certificate can be found for a given system no matter how simple its dynamics (even when a Lyapunov function is already known).

Recently, we presented a connection between the viability kernel and efficiently-computable classes of reachability constructs known as maximal reachable sets. Owing to this connection, scalable numerical algorithms (collectively referred to as Lagrangian methods [15]) such as [27, 28, 29, 30, 31, 32, 33], originally developed for maximal reachability, can now be used to approximate the viability kernel. We presented two algorithms for LTI systems with convex constraints based on piecewise ellipsoidal representations [5] and support vectors [34] that have polynomial complexity. In contrast to these results, the technique presented here reduces the complexity indirectly by decentralizing computations. The benefit of this approach is that it allows useful features of Eulerian methods such gradient-based control synthesis and handling of arbitrarily shaped nonconvex constraints be taken advantage of.

2 Problem Statement

Consider the continuous-time system


with state space (a finite-dimensional vector space), state vector , and input where is a compact (closed and bounded) and convex subset of . The vector field is assumed to be Lipschitz in and continuous in . Let


With an arbitrary, finite time horizon , for every , , and , there exists a unique trajectory that satisfies (1) and the initial condition .

For a nonempty, simply-connected, compact state constraint set we are concerned with computing the following backward construct:222By duality the arguments presented in this paper also hold for the minimal reachable tube of ; cf. [16].

Definition 1 (Viability Kernel).

The finite-horizon viability kernel333The infinite-horizon viability kernel is also known as the maximal controlled-invariant set [35]. of is the set of initial states for which there exists an input such that the trajectories emanating from those states remain in for all time :

Initial states belonging to this set are viable under (1), and the corresponding control laws are safety-preserving. The powerful Eulerian methods are capable of directly computing the viability kernel and its safety-preserving control policies. However, they rely on gridding the state space, and therefore are 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 4D in practice) due to their exponential complexity.

We restrict ourselves to LTI systems of the form


described by the matrix notation


with constant, appropriately sized and matrices.

Problem 1 (Decentralized Viability).

i) Identify a structure on and for which the viability kernel can be conservatively reconstructed from its subsystem analyses. ii) Find an isomorphic state space for (3) in which the system has this desired structure.

2.1 Preliminaries


For a set , and denote the complement and the power set of in , respectively. For brevity, denotes the infinity norm. For a constant matrix the induced norm is . For a Lebesgue measurable function defined over an interval we denote . A linear transformation of in (4) using a nonsingular matrix is defined as . A linear transformation of a set under the same mapping is .

Definition 2 (Disjoint Input).

The input is disjoint across two subsystems


of an LTI system with and if , ,


and , where is any (possibly degenerate) subset of from which the portion of the vector acting directly on subsystem draws its values.

Definition 3 (Unidirectionally Coupled).

The subsystems


with disjoint input across them are said to be unidirectionally coupled since the trajectories of (7b) are affected by those of (7a), while (7a) evolves independently from (7b). The worst-case unidirectional coupling can be characterized by .

Definition 4 (Etuc).

A subsystem is said to be externally trivially uncontrollable (ETUC) if it possesses a null input matrix.

Remark 1.

The condition on in Definition 2 enures that the inputs acting on each subsystems are independent of one another. This condition is satisfied for most physical systems where actuators are commonly uncorrelated, or for a system with an ETUC subsystem (in which case the shape of becomes irrelevant). In the most general case, however, can be (under-)approximated by a cross-product set.

3 Decentralized Viability Computation

We begin by arriving at the desired structure on system matrices that would allow for decentralized (and conservative) computation of the viability kernel. Throughout the paper we assume a partitioning of (4) that results in two subsystems. The arguments can be easily generalized to multiple subsystems as discussed in Section 4.2.

3.1 Why Decoupling of Alone is Insufficient

Consider the following system with block diagonal -matrix, and a -matrix of generic form:


Denote the two subspaces of in which the subsystems evolve as


Let be the projection of the vector onto :


and the projection of the set onto :

Lemma 1.

For any and the projection of trajectory of system (8) with initial condition is a subsystem trajectory initiating from the projection of :


. ∎

Corollary 1.

Later we will show and utilize the fact that under certain conditions this implication is bidirectional.

Proposition 1 (Wrong Approximation).

For dynamics (8) the cross-product of subsystem viability kernels of projections of is a superset of the viability kernel of :


The following counter example demonstrates that an inclusion in the opposite direction does not hold for system (8); That is, . Consider the point and constraint set . We seek to compute the viability kernel of this set under the dynamics and and input constraint . The point belongs to the cross-product of subsystem viability kernels (since subsystem 1 can use while subsystem 2 can use at the same point to keep in ), but does not belong to the actual full-order kernel (since no input exists that can keep the system in ). As such, when the system is in the form of (8) performing the analysis on subsystems would yield an over-approximation of the viability kernel. This stems from the fact that the input is non-disjoint across the subsystems. On the other hand, we do have the following correct inclusion even with a non-disjoint input.

Lemma 2.

The following holds for system (8):


Definition 5 (Ill-Posedness).

We say that a viability problem is ill-posed if the state constraint is empty.

Proposition 2 (Ill-Posed Approximation).

When is a bounded subset of (which is the case in most envelope protection problems) the approximation in Lemma 2 is ill-posed.

The proof should be clear from the fact that for any bounded set we have .

3.2 Suitable Structures for Decomposition

Consider a system with block-diagonal -matrix and a -matrix that ensures a disjoint input across the subsystems, for instance


when .

Assumption 1.

The set is a cross-product of two (arbitrarily-shaped) sets in .

Corollary 2.

Under Assumption 1 the projection of a trajectory is contained in a set if and only if the subsystem trajectories are contained in the projection of the set:

Theorem 1.

The viability kernel of under (16) can be computed exactly using subsystem kernels:

(by Assumption 1)
(via disjoint input)

Remark 2.

The use of any decomposition technique for correct (conservative) approximation of the viability kernel is contingent on satisfaction of Assumption 1 as shown previously. When does not satisfy this assumption, it can be under-approximated by the union of direct-product sets. The viability kernel can be computed for each set separately in lower dimensions (which increases the computational complexity only linearly in the number of sets). The union of the resulting kernels in full dimensions under-approximates the true viability kernel. Parallelization of viability calculations in each subspace could further reduce the computational time.

In general, we may not be able to simultaneously obtain a decoupled -matrix and a disjoint input. Instead, suppose that the system is of the form


which automatically ensures that the input is disjoint across the subsystems regardless of the shape of since one of the two (unidirectionally coupled) subsystems is ETUC (Remark 1). This system can be rewritten as


The evolution of is completely independent of the evolution of . Its effect on the lower subsystem, mapped through , can be viewed as an exogenous input to the lower subsystem, that takes values on the (possibly time-varying) subset of the upper subspace . Treating this additional input in the worst-case fashion results in conservatism. Hence, define the following construct:

Definition 6 (Discriminating Kernel).

Consider a system with adversarial inputs: control and disturbance , where is a point-wise convex and compact set-valued map from to . Let

To be conservative, we assume non-anticipative strategies for one of the inputs.444A map is non-anticipative for if for every , implies a.e.  [36]. Note that for linear systems the Isaac’s condition holds [14], and therefore it does not matter which input is selected to play with non-anticipative policies. The finite-horizon discriminating kernel of is the set of initial states for which there exists a control such that the trajectories emanating from those states remain in for every disturbance for all time :

We will use a “” subscript to distinguish a construct formed under (20) when for the lower subsystem is treated as an adversarial disturbance.

Lemma 3.

The viability kernel of a set under (20) is a superset of the discriminating kernel of when is treated as a worst-case disturbance (assumed to draw values from some time-varying set point-wise convex and compact in ) to the lower subsystem:


Let denote the trajectory of the system when is treated as a disturbance to the lower subsystem.

(a specific disturbance)

Definition 7 (Invariance Kernel).

Consider a system with a disturbance input as its only input, where is defined as in Definition 6. The finite-horizon invariance kernel of a set is the set of initial states that remain in for every disturbance for all time :

Theorem 2 (Main Decentralization Result).

The viability kernel of a set under (19) can be conservatively approximated using the subsystem viability/invariance kernels as


We first show that the inclusion holds for any set in which takes value. Since both inputs (control and “disturbance” ) are disjoint across the two subsystems we have


With , inclusion (22) follows from Lemma 3:

Note that the set-valued map at time is the finite-horizon viability kernel of the upper subsystem over the interval . This map is continuous (it is both lower and upper semicontinuous (cf. [12]) at every point in its domain) and non-decreasing [19] (i.e.  , ), with being its upper-limit in the sense of Kuratowski (Definition 8) as and its lower-limit as . Furthermore, since and are convex and compact and the dynamics linear, the sets are also convex and compact at every . From this we have that is continuous, convex and compact for every , and non-decreasing over [12].

We use these statements to argue that a digression from the formulation in (22) loses its sufficiency to guarantee an under-approximation in the sense that if the uncertainty set is assumed to be a subset of for any then the cross-product may not generate an under-approximation of the viability kernel: Consider a set-valued map s.t. , (e.g. a constant set ). It is clear from (23) that


since for any set , and therefore . There is no guarantee that this superset in (24) is a subset of ; Lemma 3 is no longer applicable. On the flip side, if is such that for any (e.g. a constant set ), then an excessively conservative under-approximation is obtained. ∎

3.3 Sub-Interval Formulation and Decentralized Algorithm

In practice, we can perform the analysis over sub-intervals (similarly to [37]) while still maintaining conservatism. During each sub-interval the set is sampled and kept constant in backward time. Such sub-interval analysis is possible via the semi-group property in both subspaces as well as the following results in .

Proposition 3.

For , time steps each of length we have that


where with .


Notice that since is a non-decreasing sequence of compact and convex sets with we have that for a fixed , for every , . Using this, the fact that , and the semi-group property we have

where is the concatenation of functions over .

In the limit this set converges to the invariance kernel with unsampled input set.

Definition 8 (Kuratowski upper and lower limits [19]).

Let be a sequence of subsets in a metric space . The upper-limit of as is

where . Its lower-limit is

Proposition 4.

Denote by the intersection of sub-interval invariance kernels from Proposition 3. For the sequence of subsets we have


Given , define a piecewise constant set-valued map for which is the unique integer in satisfying when varies backwards from to (i.e. a backward sample and hold of ). Recall that is non-decreasing and continuous, and compact for every . Clearly, . The sequence converges to from outside: We say that iff . As , , where denotes the ball (associated with a metric ) of radius centered at . In other words, , s.t. . So , and therefore . On the other hand, we know from the semi-group property that . Hence,

Using this formulation we can perform the decentralized analysis in Theorem 2 via Algorithm 1 over sub-intervals.

1: Assumed integer.
4:for  to  do
7:end for
Algorithm 1 Sub-Interval Decentralized Computations

3.4 Bounding the Approximation in

Notice from Theorem 2 that the computed construct in the upper subspace is exact in that


On the other hand additional conservatism is introduced in the lower subspace due to treating the effect of the upper subsystem as a worst-case disturbance. Quantifying this error remains an open problem. However, we can formulate a qualitative lower bound on the shrinkage of the invariance kernel in in backward time. This bound will be expressed in terms of system-specific (and ultimately, design-specific) parameters that form the desired structure (19):

Following [38], the invariance kernel in can be expressed as


with denoting the Pontryagin difference. Let be the norm-ball of radius about the origin, and define ,


Bounding the contribution of the uncertainty (disturbance) in computation of the invariance kernel over the interval we have [37] that


Clearly, this contribution is weakened as . Further, we have


with . From the dual of the results in [37], we know that the Hausdorff distance of the two sets in the inclusion above decreases as , and tends to zero if . The Kuratowski upper-limit of the left-hand-side of (34) is therefore as (via Proposition 4). Now, notice that for sufficiently small ,


where and respectively denote the largest singular value and the dimension of the lower subsystem. Therefore (34) provides a qualitative lower-bound on how much can shrink in backward time in terms of , the magnitude of the unidirectional coupling , the supremum of (the viability kernel in ), and the largest singular value of the lower subsystem. If we can choose