Linear Reachability

Chapter 1: Forward Reachability

A system starts in some initial state x₀ drawn from an initial set X₀. At each time step, the state evolves according to the dynamics, subject to disturbances from a disturbance set W. The reachable set at time t is the set of all states the system could possibly be in, over all possible initial conditions and all possible disturbance sequences:

This is a set, not a distribution. It contains every state that could result from any combination of allowed initial conditions and disturbances. If R_t ∩ X_unsafe = ∅ for all t, the system is safe.

The key recursion for computing the reachable set one step at a time is:

Start with R₀ = X₀. At each step, apply the dynamics f to every state in R_t and every disturbance in W. The result is R_t+1. This is conceptually simple but computationally brutal for general f — we'd need to track an infinite set of points. The power of linear dynamics is that this infinite operation reduces to just two set operations.

Chapter 2: Linear Systems

For a general nonlinear system x_t+1 = f(x_t, w_t), computing the reachable set is intractable — it requires tracking the image of a set under an arbitrary function. But for linear time-invariant (LTI) systems, the dynamics have the form:

where A is the state transition matrix and B maps disturbances into state space. Both A and B are constant matrices. The initial set X₀ and disturbance set W are compact convex sets (think: ellipsoids, boxes, or polytopes).

Substituting into the reachability recursion:

This decomposes into two operations:

The symbol ⊕ denotes the Minkowski sum: P ⊕ Q = {p + q : p ∈ P, q ∈ Q}. We will define both operations precisely in the next chapter. The crucial point is that these two operations — linear transformation and Minkowski sum — are the only tools we need. Everything else in this chapter is about computing them efficiently.

Let's trace one step in 2D. Suppose R₀ is a unit square centered at the origin, A = [[0.9, 0.1], [−0.1, 0.9]] (a slight rotation + contraction), and BW is a small box [−0.1, 0.1]².

Each step rotates and shrinks (because A has eigenvalues with magnitude 0.9), then fattens by the disturbance. Over many steps, the reachable set converges to a limit set — the disturbance injection balances the contraction. If that limit set avoids the unsafe region, the system is safe forever.

Chapter 3: Set Operations

Let's define the two building blocks precisely and build intuition for each.

Linear transformation. Given a set P and a matrix A, the image is:

Geometrically, A stretches, rotates, and/or shears every point in P. If P is a circle, AP is an ellipse. If P is a square, AP is a parallelogram. The key property: if P is convex, AP is convex.

Minkowski sum. Given two sets P and Q:

Geometrically: slide Q around the boundary of P (or equivalently, slide P around Q). The swept area is P ⊕ Q. Equivalently, "fatten" P by the shape of Q. If P is a square and Q is a disk, P ⊕ Q is a rounded rectangle.

The computational cost of these operations depends heavily on how we represent the sets. A set described by 4 vertices is easy to transform; a set described by 1000 vertices is harder. This brings us to the critical question: what data structure should we use for sets? The next three chapters explore three answers — polytopes, zonotopes, and support functions — each with different trade-offs.

Chapter 4: The Dependency Effect

Before we choose a set representation, we need to understand a subtle trap. Set operations treat their operands as independent. Consider computing A · (P ⊕ Q). If we first compute S = P ⊕ Q, then compute A · S, the result is correct. But what about A · P ⊕ A · Q? Is that the same?

For linear maps, yes: A(P ⊕ Q) = AP ⊕ AQ. This is the distributive property. But trouble arises when the same variable appears multiple times in an expression. Consider:

But if we treat the two x's independently:

In 1D it happens to agree. But in higher dimensions, it often doesn't. Consider x ∈ [−1,1]² and the expression Ax + Bx with A and B being different matrices. The true result is {(A+B)x : x ∈ [−1,1]²}, but Minkowski-summing AP ⊕ BP treats the x in each term independently, potentially producing a larger set. This is the dependency effect (also called the wrapping effect).

The practical implication: always simplify expressions before applying set operations. Compute (A+B)x as a single linear map on x, rather than computing Ax and Bx separately and summing. This algebraic preprocessing is critical for controlling wrapping.

In the reachability recursion R_t+1 = A · R_t ⊕ BW, the dependency effect is manageable because x_t and w_t really are independent (the current state and the current disturbance are separate variables). The danger arises when representing R_t itself involves implicit variables that get duplicated across operations.

Chapter 5: Polytopes

A polytope is the most natural way to represent a convex set: it's a bounded region defined by flat faces. In 2D, polytopes are convex polygons. In 3D, convex polyhedra. There are two equivalent representations.

H-polytope (half-space representation): the intersection of finitely many half-spaces.

Each row of H defines a face normal, and the corresponding entry of h defines the face offset. A cube in 3D needs 6 half-spaces. An n-gon in 2D needs n half-spaces.

V-polytope (vertex representation): the convex hull of finitely many points.

The set of all convex combinations of the vertices. A triangle has 3 vertices. A cube has 8.

Each representation has a "cheap" operation and an "expensive" one:

The critical problem for reachability is the Minkowski sum. In V-representation, if P has m vertices and Q has n vertices, P ⊕ Q can have up to O(mn) vertices in 2D and exponentially more in higher dimensions. Since reachability requires a Minkowski sum at every time step, the vertex count grows exponentially with the number of steps.

Chapter 6: Zonotopes

A zonotope is a special type of polytope defined by a center and a set of generators. Given a center c ∈ R^d and generator matrix G = [g₁, …, g_k] (each g_i ∈ R^d), the zonotope is:

Each generator "stretches" a unit interval along its direction. The zonotope is the Minkowski sum of all these line segments, translated by c. In 2D, a zonotope with 2 generators is a parallelogram. With 3 generators, it's a hexagon. With k generators, it's a convex polygon with at most 2k sides.

Let's trace the reachability step. Suppose R_t = Z(c_t, G_t) with k_t generators, and BW = Z(0, G_w) with k_w generators.

After one step: k_t+1 = k_t + k_w generators. After T steps: k₀ + T · k_w generators. Linear growth! If the initial set has 2 generators and the disturbance has 2 generators, after 100 steps we have 202 generators — not 10³⁰ vertices.

Chapter 7: Overapproximation

Exact reachability is ideal but often too expensive. Overapproximation computes a set P̃ that contains the true reachable set: R_t ⊆ P̃_t. If the overapproximation doesn't intersect the unsafe region, neither does the true reachable set. Safety is preserved.

Common overapproximation strategies:

The critical insight: overapproximation introduces one-sided error. If P̃ ∩ X_unsafe ≠ ∅, we cannot conclude the system is unsafe — the intersection might be entirely within the overapproximation gap. We can only conclude safety (when there's no intersection) or "inconclusive" (when there is). To resolve inconclusive cases, we refine the overapproximation or switch to a tighter method.

Chapter 8: Support Functions

All the set representations so far store explicit geometric data (vertices, generators, half-spaces). The support function takes a fundamentally different approach: it represents a convex set by answering the question "how far does the set extend in direction d?" for any direction d.

The support function ρ_P(d) gives the maximum "projection" of P along direction d. It completely characterizes any convex compact set — knowing ρ_P(d) for all directions d uniquely determines P.

Support functions compose beautifully with our two operations:

The Minkowski sum becomes simple addition! No vertex explosion, no generator concatenation. The cost is that the support function is defined for all directions d, not just a finite set. In practice, we evaluate it at a finite set of directions d₁, …, d_n and construct an outer-approximating polytope:

More directions = tighter approximation. The directions can be chosen adaptively — place more directions where the boundary has high curvature.

Chapter 9: LP Safety Checking

We now have tools to compute (or overapproximate) the reachable set. The final question is: does it intersect the unsafe region? This is a safety checking problem, and for polytopic sets, it reduces to linear programming (LP).

Suppose the unsafe set is a polytope X_unsafe = {x : Cx ≤ d}. And our (over-approximated) reachable set P̃ is also a polytope or can be queried via a support function. Safety holds if and only if P̃ ∩ X_unsafe = ∅.

To check intersection, we can test each face of X_unsafe. Face j is defined by c_j^Tx ≤ d_j. The reachable set violates this constraint if and only if:

If this holds for at least one face j — meaning the reachable set extends beyond at least one face of the unsafe set in the outward direction — then P̃ might not actually intersect X_unsafe. Wait, that's the wrong direction. Let's be precise:

Alternatively, using support functions: P̃ and X_unsafe are disjoint if and only if there exists a separating hyperplane. By the separating hyperplane theorem, this is equivalent to finding a direction d such that:

In practice, the LP approach is most common. Modern LP solvers handle millions of constraints in seconds, making safety checking fast even for complex reachable sets.

pseudocode
# Safety verification via reachability
Given: A, B, X0, W, Xunsafe, time horizon T

R ← X0
for t = 1 to T:
  R ← overapprox(A · R ⊕ B · W)   # zonotope/support fn step
  if LP_feasible(R ∩ Xunsafe):
    return "INCONCLUSIVE at step t"
return "SAFE for [0, T]"

The full pipeline: choose a set representation (zonotope for speed, support function for tightness), propagate the reachable set through the linear dynamics, overapproximate at each step to control complexity, and check intersection with the unsafe set via LP. If the LP is infeasible at every step, the system is provably safe. If feasible at any step, refine the overapproximation or increase the number of support directions.