Ellipsoidal Techniques for Reachability Analysis of ... - Semantic Scholar
Ellipsoidal Techniques for Reachability Analysis of Discrete-Time Linear Systems∗ Alex Kurzhanskiy and Pravin Varaiya Department of Electrical Engineering and Computer Sciences University of California, Berkeley, CA June 30, 2005
Abstract This paper describes the computation of reach sets for discrete-time linear control systems with time-varying coefficients and ellipsoidal bounds on the controls and initial conditions. The algorithms construct external and internal ellipsoidal approximations so that they touch the reach set boundary from outside and from inside. Recurrence relations that describe the time evolution of these approximations are provided. The paper also deals with discrete-time linear systems with singular state transition matrix.
1
Introduction
We extend the use of ellipsoidal methods for continuous-time linear systems developed in [1] to calculate the reach set for discrete-time linear systems. We also deal with systems with a singular state transition matrix and show how to approximate the reach set with any given accuracy using regularization and provide the means for picking the appropriate regularization parameters that guarantee this accuracy over a specified finite time interval. We make no controllability assumptions regarding the system. Thus, the suggested regularization technique also extends results of [1] to systems that are not completely controllable. In case of singular state transition matrix, or if the system is not controllable at every time step, the exact representation of the reach set on every time step by external ellipsoidal approximation is not possible. But the exact representation is possible for the regularized reach set that overapproximates the actual one, and whose boundary can be made arbitrarily close to the boundary of the actual one over any given finite time interval. Using internal ellipsoidal approximations, on the other hand, it is always possible to exactly represent ∗
Research supported by NSF Grant CCR-00225610.
1
the actual reach set. However, the concept of ‘good curves’ that was introduced in [1] and without which the computation process becomes heavy and complicated, can only be applied for systems with nonsingular state transition matrix. Therefore, in order to make use of good curves in the singular case, we perform the internal approximation of the regularized reach set.
2
Related work
Research on dynamical and hybrid systems has produced methods for verification and controller synthesis. A common step in most methods is the calculation of reach sets. For example, [2] uses bounded model checking for verification of discrete-time linear hybrid systems. For that, the hybrid system is converted into a discrete transition system. The conversion requires the computation of reach sets of discrete-time linear systems. Exact computation of reach sets is not generally possible, and approximation techniques are needed. The choice of the reach set representation determines the efficiency of such techniques. The more complex the representation is, the more costly is the storage of the sets, the more elaborate is the computation, but the better is the reach set approximation. Choosing the reach set representation is a compromise between these factors. The level set method [3, 4] deals with general nonlinear controlled systems and gives exact representation of their reach sets, but the computation is complex, requiring solving the HJB equation and finding the set of states that belong to sub-zero level set of the value function. The method [4] is impractical for systems of dimension higher than three. Reachability analysis of controlled linear systems with convex bounds on the control and initial conditions reduces to effective manipulation with convex sets, performing set-valued operations such as unions, intersections, geometric sums and differences. Two basic objects are used as convex approximations: various kinds of polytopes, e.g. general polytopes, zonotopes, parallelotopes, rectangular polytopes; and ellipsoids. Reachability analysis and verification for discrete-time linear systems through general polytopes are implemented in the Multi Parametric Toolbox (MPT) for Matlab [5], [6]. The reach set at every time step is computed as the geometric sum of two polytopes. The procedure consists in finding the vertices of the resulting polytope and calculating its convex hull. The convex hull algorithm employed by MPT is based on the Double Description method [7] and implemented in the CDD/CDD+ package [8]. Its complexity is V n , where V is the number of vertices and n is the state space dimension. Hence the use of MPT is practical for low dimensional systems. But even in low dimensional systems the number of vertices in the the reach set polytope can grow very large with the number of time steps. For example, consider the system, xk+1 = Axk + uk ,
2
cos 1 − sin 1 with A = , uk ∈ {u ∈ R2 | u∞ ≤ 1}, and x0 ∈ {x ∈ R2 | x∞ ≤ 1}. sin 1 cos 1 Starting with a rectangular initial set, the number of vertices of the reach set polytope is 4k + 4 at the kth step. In d/dt [9], the reach set is approximated by unions of rectangular polytopes [10]. The x 1(t 2)
R (t 2) x 2(t 2) R [t 1,t 2] R (t 1) x 1(t 1)
x 2(t 1)
x 1(t 1)
x 2(t 1)
R [t 0,t 1] x 1(t 0)
x 2(t 0)
(a)
x 1(t 0)
x 2(t 0)
(b)
(c)
(d)
(e)
(f)
Figure 1: Reach set approximation by union of rectangles. algorithm works as follows. First, given the set of initial conditions defined as a polytope, the evolution in time of the polytope’s extreme points is computed (figure 1(a)). R(t1 ) in figure 1(a) is the reach set of the system at time t1 , and R[t0 , t1 ] is the set of all points that can be reached during [t0 , t1 ]. Second, the algorithm computes the convex hull of vertices of both, the initial polytope and R(t1 ) (figure 1(b)). The resulting polytope is then bloated to include all the reachable states in [t0 , t1 ] (figure 1(c)). Finally, this overapproximating polytope is in its turn overapproximated by the union of rectangles (figure 1(d)). The same procedure is repeated for the next time interval [t1 , t2 ], and the union of both rectangular approximations is taken (figure 1(e,f)), and so on. Rectangular polytopes are easy to represent and the number of facets grows linearly with dimension, but a large number of rectangles must be used to assure the approximation is not overly conservative. Besides, the important part of this method is again the convex hull calculation whose implementation relies on the same CDD/CDD+ library. This limits the dimension of the system and time interval for which it is practical to calculate the reach set. Polytopes can give arbitrarily close approximations to convex sets, but the number of vertices can grow prohibitively large and, as shown in [11], the computation of a polytope by its convex hull becomes intractable for large number of vertices in high dimensions. The method of zonotopes for external approximation of reach sets is presented in [12]. A 3
zonotope is a special class of polytopes (see [13]) of the form, n
Z = {x ∈ R | x = c +
p
αi gi , − 1 ≤ αi ≤ 1},
i=1
wherein c and g1 , ..., gp are vectors in Rn . Thus, a zonotope Z is defined by its center c and generator vectors g1 , ..., gp . The value p/n is called the order of the zonotope. This compact representation and the fact that they are closed under linear transformation and geometric sum, make zonotopes attractive for reach set computation. The difficulty is that with every integration step the order of the approximating zonotope increases by p/n. The difficulty can be averted by limiting the number of generator vectors, and overapproximating zonotopes whose number of generator vectors exceeds this limit by lower order zonotopes. The zonotope method of [12] provides neither a notion of tightness of approximation nor an estimate of the error resulting from the zonotope order reduction. Further, the method does not indicate any way of control synthesis, i.e. finding a sequence of controls that drives the system to the desired terminal set or point. CheckMate [14] is a Matlab toolbox that can evaluate specifications for trajectories starting from the set of initial (continuous) states corresponding to the parameter values at the vertices of the parameter set. This provides preliminary insight into whether the specifications will be true for all parameter values. The method of oriented rectangluar polytopes for external approximation of reach sets is introduced in [15]. The basic idea is to construct an oriented rectangular hull of the reach set for every time step, whose orientation is determined by the singular value decomposition of the sample covariance matrix for the states reachable from the vertices of the initial polytope. The limitation of CheckMate and the method of oriented rectangles is that only autonomous (i.e. uncontrolled) systems are allowed, and only an external approximation of the reach set is provided. Requiem [16] is a Mathematica notebook which, given a linear system, the set of initial conditions and control bounds, symbolically computes the reach set exactly, using the experimental quantifier elimination package. Quantifier elimination is the removal of all quantifiers (the universal quantifier ∀ and the existential quantifier ∃) from a quantified system. Each quantified formula is substituted with quantifier-free expression with operations +, ×, = and k0 from the initial position (k0 , x0 ) is the set of all states x[k] reachable at time k by the system (1) with x[k0 ] = x0 through all possible controls that satisfy (3). The reach set X (k, k0 , X0 ) is X (k, k0 , X0 ) =
{X (k, k0 , x0 ) | x0 ∈ X0 )}.
Some important facts about reach sets follow. Lemma 3.1 The set-valued map X (k, k0 , X0 ) satisfies the semigroup property X (k, k0 , X0 ) = X (k, i, X (i, k0 , X0 )),
k0 ≤ i ≤ k.
(5)
Lemma 3.2 The reach set X (k, k0 , E(x0 , X0 )) can be expressed as X (k, k0 , E(x0 , X0 )) = q[k] + Φ(k, k0 )E(0, X0 ) +
k−1
Φ(k, i + 1)E(0, B[i]P [i]B T [i]),
(6)
i=k0
where q[k] = Φ(k, k0 )x0 +
k−1
Φ(k, i + 1)B[i]p[i].
(7)
i=k0
Using (6) and (4), we obtain the next result. Lemma 3.3 The support function of the reach set is ρ(l | X (k, k0 , E(x0 , X0 ))) = l, q[k] + l, Φ(k, k0 )X0 ΦT (k, k0 )l 1/2 +
k−1
l, Φ(k, i + 1)B[i]P [i]B T [i]ΦT (k, i + 1)l 1/2 .
i=k0
Representation (6) leads to the next fact. 6
(8)
Lemma 3.4 The reach set X (k, k0 , E(x0 , X0 )) is a convex compact set in Rn . If the matrices A[i], k0 ≤ i < k, are nonsingular, the ellipsoid Φ(k, k0 )E(0, X0 ) in sum (6) is nondegenerate, and so the reach set X (k, k0 , E(x0 , X0 )) contains an open set of Rn . Its boundary points have an important characterization: y is a boundary point of X (k, k0 , E(x0 , X0 )) only if there exists a support vector ly = 0 so that ly , y = ρ(ly | X (k, k0 , E(x0 , X0 ))) = max{ ly , x | x ∈ X (k, k0 , E(x0 , X0 ))}.
(9)
The control uy [i], k0 ≤ i < k, and the initial state x[k0 ] = x0y ∈ E(x0 , X0 ), which drive the system (1) from the state x[k0 ] = x0y to a boundary point x[k] = y is determined from the following theorem. Theorem 3.1 Suppose A[i], k0 ≤ i < k, are nonsingular, and x[k] = y is a boundary point of X (k, k0 , E(x0 , X0 )). Then the control uy [i], k0 ≤ i < k, and the initial state x[k0 ] = x0y , which yield the unique trajectory xy [i] that reaches y = xy [k] from x0y = x[k0 ], and ensures (9), satisfy the following maximum principle for control: ly , Φ(k, i + 1)B[i]uy [i] = max{ ly , Φ(k, i + 1)B[i]u | u ∈ Φ(k, i + 1)E(p[i], P [i])} = ly , Φ(k, i + 1)B[i]p[i]
+ ly , Φ(k, i + 1)B[i]P [i]B T [i]ΦT (k, i + 1)ly 1/2
(10)
for k0 ≤ i < k, and the maximum condition for the initial state: ly , Φ(k, k0 )x0y = max{ ly , x | x ∈ Φ(k, k0 )E(x0 , X0 )} = ly , Φ(k, k0 )x0 + ly , Φ(k, k0 )X0 ΦT (k, k0 )ly 1/2 .
(11)
Here ly is the support vector for the set X (k, k0 , E(x0 , X0 )) at the point y that satisfies (9). Remark. To simplify the notation we shall denote R[i] = B[i]P [i]B T [i] and from now on use R[i] instead of B[i]P [i]B T [i] where possible.
4
External ellipsoidal approximations
Relation (6) shows that the reach set X (k, k0 , E(x0 , X0 )) is a sum of k − k0 + 1 ellipsoids. We overapproximate the reach set by a single ellipsoid. Here we consider the situation when matrices A[i] are nonsingular. The case with singular A[i] is treated in section 6. In the beginning we shall also assume that matrices R[i] are nonsingular. Remark. R[i] are nonsingular if the matrices B[i] ∈ Rn×m have rank n (n ≤ m). We shall look for external approximating ellipsoids for the reach set X (k, k0 , E(x0 , X0 )) in the class of ellipsoids E(q[k], Q[k]) that satisfy the following two relations: q[k] = A[k − 1]q[k − 1] + B[k − 1]p[k − 1], 7
q[k0 ] = x0 ,
(12)
and Q[k] =
k−1
si [k] + s0 [k]
i=k0
×
k−1 i=k0
T s−1 i [k]Φ(k, i + 1)R[i]Φ (k, i + 1)
T +s−1 0 [k]Φ(k, k0 )X0 Φ (k, k0 ) ,
Q[k0 ] = X0 ,
(13)
for arbitrary s0 [k] > 0, si [k] > 0 for k > k0 and k0 ≤ i < k. Using (8), it can be easily checked that ρ(l | E(q[k], Q[k])) ≥ ρ(l | X (k, k0 , E(x0 , X0 ))) for all directions l. Hence, X (k, k0 , E(x0 , X0 )) ⊆ E(q[k], Q[k]). Definition 4.1 We say that ellipsoid E(q[k], Q[k]), defined by (12) and (13), is a tight external approximation of the reach set X (k, k0 , E(x0 , X0 )) if there exists direction l such that ρ(±l | E(q[k], Q[k])) = ρ(±l | X (k, k0 , E(x0 , X0 ))). A tight external ellipsoid defined by direction l, touches the boundary of the reach set at point xl such that ρ(l | X (k, k0 , E(x0 , X0 ))) = l, xl = ρ(l | E(q[k], Q[k])) = l, q[k] + l, Q[k]l 1/2 . This condition will be fulfilled if we choose
and
si [k] = l, Φ(k, i + 1)R[i]ΦT (k, i + 1)l 1/2 ,
(14)
s0 [k] = l, Φ(k, k0 )X0 ΦT (k, k0 )l 1/2 .
(15)
To check it, we substitute (14) and (15) into (13), calculate the support function of the resulting ellipsoid and arrive at (8). As we can see, for every direction l there exists a tight external ellipsoidal approximation of the reach set X (k, k0 , E(x0 , X0 )), whose shape matrix is calculated by (13). For all 0 ≤ i < k, the parameters si [k] in (15) and (14) depend on k. As a result, the sum (13) needs to be recomputed for every k with new values of s0 [k] and si [k]. To get rid of this dependency on k, and thus greatly reduce the computation effort, we introduce the notion of ‘good curves’. These are trajectories l[k] that satisfy the homogeneous adjoint equation l[k] = AT [k]l[k + 1], 8
l(k0 ) = l0 .
(16)
Vector l[k] can be expressed as l[k] = ΦT (k0 , k)l0
(for time-invariant systems: l[k] = ((Ak−k0 )T )−1 l0 ).
(17)
With l[k] defined in (17), relations (14) and (15) simplify into si = l0 , Φ(k0 , i + 1)R[i]Φ(k0 , i + 1)l0 1/2
(18)
s0 = l0 , X0 l0 1/2 .
(19)
and
Now si and s0 no longer depend on k. When substituted into (13), they give us a recurrence relation for matrix Q[k], Q[k + 1] = (1 + π[k])A[k]Q[k]AT [k] + (1 + π −1 [k])R[k],
Q[k0 ] = X0 ,
(20)
in which π[k] is defined by sk i=k0 si + s0
π[k] =
k−1
=
k−1
=
l0 , Φ(k0 , k + 1)R[k]ΦT (k0 , k + 1)l0 1/2 . l0 , Φ(k0 , k)Q[k]ΦT (k0 , k)l0 1/2
l0 , Φ(k0 , k + 1)R[k]ΦT (k0 , k + 1)l0 1/2 T 1/2 + l , X l 1/2 0 0 0 i=k0 l0 , Φ(k0 , i + 1)R[i]Φ (k0 , i + 1)l0
(21)
The expression for the center of the approximating ellipsoid is given by q[k + 1] = A[k]q[k] + B[k]p[k],
q[k0 ] = x0 .
(22)
In both (20) and (22) k ≥ k0 . We state this result as a theorem. Theorem 4.1 Let l[k] satisfy the equation (17). Then the ellipsoid E(q[k], Q[k]) given by (20) and (22), is such that for all k ≥ k0 the reach set X (k, k0 , E(x0 , X0 )) ⊆ E(q[k], Q[k]) and ρ(l[k] | E(q[k], Q[k])) − ρ(l[k] | X (k, k0 , E(x0 , X0 ))) = 0. The supporting hyperplane for the reach set X (k, k0 , E(x0 , X0 )) generated by the vector l[k] is also a supporting hyperplane for E(q[k], Q[k]) and touches it at the same point xl [k], given by Q[k]l[k] . (23) xl [k] = q[k] + l[k], Q[k]l[k] 1/2 Since for every boundary point xl [k] of the reach set X (k, k0 , E(x0 , X0 )) there exists direction l[k] and corresponding tight external ellipsoid E(q[k], Q[k]), the following result is established. 9
Theorem 4.2 The reach set X (k, k0 , E(x0 , X0 )) can be expressed as X (k, k0 , E(x0 , X0 )) =
E(q[k], Q[k]).
(24)
l[k]
If l[k] is a “good curve”, it specifies the external ellipsoid E(q[k], Q[k]), with Q[k] and q[k] coming from (20) and (22), which touches the reach set at some point xl [k]. In order for the points xl [·] to form the system trajectory on k0 , ..., k, that is, for xl [k] to satisfy xl [k] = Φ(k, k0 )x0 +
k−1
Φ(k, i + 1)B[i]u[i],
i=k0
the corresponding control u[i] ∈ E(p[i], P [i]) and the initial condition x0 ∈ E(x0 , X0 ) are derived from the maximum principle for control (10) and maximum condition (11). They are specified by P [i]B T [i]ΦT (k0 , i + 1)l0 + p[i], (25) u[i] = l0 , Φ(k0 , i + 1)R[i]ΦT (k0 , i + 1)l0 1/2 and x0l =
X0 l0 + x0 . l0 , X0 l0 1/2
(26)
Thus, if the system starts at x0l at time step k0 , the control u[i] will bring it at time step k to the boundary point of the reach set, at which the reach set is touched by the external approximating ellipsoid E(q[k], Q[k]) defined by the direction l[k]. Remark. For n-dimensional system, every step of computation of external ellipsoidal approximation requires 2n2 scalar multiplications for the center of ellipsoid and (8n3 + 4n2 + 2n) scalar multiplications for its shape matrix. Since the trajectory of the center is the same for different values of parameter l0 , computing L ellipsoidal approximations for k time steps requires k(L(8n3 + 4n2 + 2n) + 2n2 ) scalar multiplications. So, the complexity of computation grows linearly with the number of time steps and number of approximations, and polynomially, as O(n3 ), with the state space dimension.
5
Singular R
We start by stating a simple fact that will be useful in our proofs. Lemma 5.1 Let a ≥ b > 0. Then from a2 − b2 < ε2 for some ε > 0, it follows that a − b < ε. 10
In case the matrix R[i] is singular for some i, the recurrence relations (20) are not valid because we cannot guarantee that parameter si defined in (18) is strictly positive. Thus, we need to replace R[i] with a suitable nonsingular matrix. Define Rα [i], Rα [i] = R[i] + α2 I,
k0 ≤ i < k,
(27)
for α > 0. Rα [i] is positive definite and Rα [i] → R[i] as α → 0. Corresponding to Rα [i] we define the “regularized” reach set Xα (k, k0 , E(x0 , X0 )), Xα (k, k0 , E(x0 , X0 )) = q[k] + Φ(k, k0 )E(0, X0 ) +
k−1
Φ(k, i + 1)E(0, Rα [i]),
(28)
i=k0
with support function ρ(l | Xα (k, k0 , E(x0 , X0 ))) = l, q[k] + l, Φ(k, k0 )X0 ΦT (k, k0 )l 1/2 +
k−1
l, Φ(k, i + 1)Rα [i]ΦT (k, i + 1)l 1/2 ,
(29)
i=k0
with q[k] determined from (7). The difference between the support functions of the regularized reach set Xα (k, k0 , E(x0 , X0 )) and the actual reach set X (k, k0 , E(x0 , X0 )) is ρ(l | Xα (k, k0 , E(x0 , X0 ))) − ρ(l | X (k, k0 , E(x0 , X0 ))) =
k−1
( l, Φ(k, i + 1)Rα [i]ΦT (k, i + 1)l 1/2 − l, Φ(k, i + 1)R[i]ΦT (k, i + 1)l 1/2 ).
i=k0
Let us estimate the ith element of this sum, l, Φ(k, i + 1)Rα [i]ΦT (k, i + 1)l 1/2 − l, Φ(k, i + 1)R[i]ΦT (k, i + 1)l 1/2 ,
k0 ≤ i < k.
Let σ(A[i]) be the largest singular value of matrix A[i], and denote σmax = max {σ(A[i])}. k0 ≤i 0 there exists α > 0, such that 0 ≤ ρ(l | Xα (k, k0 , E(x0 , X0 ))) − ρ(l | X (k, k0 , E(x0 , X0 ))) < ε for all k ≥ k0 , and all vectors l. More precisely, α can be chosen as ε , α= k−k0 −1 2 2(σmax + · · · + σmax + σmax + 1) where σmax is defined in (30). 11
(31)
Remark. In (30) we could have taken i strictly greater than k0 . We leave this definition as is, however, because it will be used in the next section in its present form. Substituting R with Rα in (18) and (21), we obtain recurrence relations for Qα [k], the shape matrix of the tight external ellipsoid of the regularized reach set Xα (k, k0 , E(x0 , X0 )). The recurrence relation for the center of this ellipsoid (22) remains untouched. So, the relation X (k, k0 , E(x0 , X0 )) =
E(q[k], Qα [k])
(32)
l[k]
holds with ε accuracy in the Hausdorff metric. Expression (25) for the control is similarly modified by replacing R with Rα .
6
Singular A
When A[i] is singular for some i, k0 ≤ i < k, relation (17) is invalid because the inverse for the state transition matrix Φ(k, k0 ) does not exist; hence, we cannot use the expression Φ(k0 , k); and the subsequent formulas (18) and (20) do not make sense for singular A[i]. The way around it is to regularize matrix A[i]. Let the singular value decomposition of A[i] be A[i] = Ui Σi Vi , and define the new matrix Aδ [i] by Aδ [i] = Ui (Σi + δI)Vi = A[i] + δUi Vi ,
k0 ≤ i < k,
(33)
where δ > 0. As one can see, Aδ [i] is nonsingular for δ > 0, and Aδ [i] → A[i] as δ → 0. Corresponding to the matrix Aδ [i] there is an invertible state transition matrix Φδ (k, k0 ). Remark. Here, actually, δ may depend on i, and Φδ (k, k0 ) is Φδ (k, k0 ) = Aδk−1 [k − 1]Aδk−2 [k − 2] · · · Aδk0 +1 [k0 + 1]Aδk0 [k0 ]. We drop the subindices, however, to simplify the notation. Let H be any symmetric positive definite matrix. We would like to find such δ > 0 that for given ε > 0, −ε2 < l, Aδ [i]HATδ [i]l − l, A[i]HAT [i]l < ε2 . for all vectors l, and any i, then by lemma 5.1, −ε < l, Aδ [i]HATδ [i]l 1/2 − l, A[i]HAT [i]l 1/2 < ε, For that we write Aδ [i]HATδ [i] − A[i]HAT [i] = δ(A[i]HViT UiT + Ui Vi HAT [i]) + δ2 Ui Vi HViT UiT . 12
If we denote the largest singular value of H as h, and the largest singular value of the symmetric matrix A[i]HViT UiT + Ui Vi HAT [i] as λ = λ(A[i]HViT UiT + Ui Vi HAT [i]), then δ could be found from the inequality hδ2 + λδ < ε2 . Thus, we arrive at the following lemma. Lemma 6.1 Let H = H T > 0. Then, for any ε > 0 there exists δ > 0 such that −ε < l, Aδ [i]HATδ [i]l 1/2 − l, A[i]HAT [i]l 1/2 < ε. More precisely, δ can be chosen as
2
−λ + λ + 4ε2 , δ= 2h + 1
(34)
where λ is the largest singular value of the matrix A[i]HViT UiT + Ui Vi HAT [i]. We define the regularized reach set Xδ,α (k, k0 , E(x0 , X0 )) by Xδ,α(k, k0 , E(x0 , X0 )) = q[k] + Φδ (k, k0 )E(0, X0 ) +
k−1
Φδ (k, i + 1)E(0, Rα [i]).
(35)
i=k0
Its support function is ρ(l | Xδ,α (k, k0 , E(x0 , X0 ))) = l, q[k] + l, Φδ (k, k0 )X0 ΦTδ (k, k0 )l 1/2 +
k−1 i=k0
l, Φδ (k, i + 1)Rα [i]ΦTδ (k, i + 1)l 1/2 ,
(36)
where Rα [i] is defined in (27) and q[k] comes from (7) as before. Our goal is to make the regularized reach set arbitrary close to the actual reach set for every k ≥ k0 , by the appropriate choice of parameters δi , k0 ≤ i < k and α. In other words, for any ε > 0 we should find δi > 0, k0 ≤ i < k, and α > 0, such that 0 ≤ ρ(l | Xδ,α (k, k0 , E(x0 , X0 ))) − ρ(l | X (k, k0 , E(x0 , X0 ))) < ε for all k ≥ k0 , and all vectors l. We now show how to construct such α and the sequence δi . Let us estimate the difference of the support functions of the regularized reach set and the actual reach set ρ(l | Xδ,α (k, k0 , E(x0 , X0 ))) − ρ(l | X (k, k0 , E(x0 , X0 ))) = ρ(l | Aδ [k − 1]Xδ,α (k − 1, k0 , E(x0 , X0 ))) − ρ(l | A[k − 1]Xδ,α (k − 1, k0 , E(x0 , X0 ))) +ρ(l | A[k − 1]Xδ,α (k − 1, k0 , E(x0 , X0 ))) − ρ(l | A[k − 1]X (k − 1, k0 , E(x0 , X0 ))) +ρ(l | E(0, Rα [k − 1])) − ρ(l | E(0, R[k − 1])). 13
(37)
As we can see, this difference breaks up into three components. First, we analyze ρ(l | Aδ [k − 1]Xδ,α (k − 1, k0 , E(x0 , X0 ))) − ρ(l | A[k − 1]Xδ,α (k − 1, k0 , E(x0 , X0 ))) =
k−2 i=k0
( l, Aδ [k − 1]Φδ (k − 1, i + 1)Rα [i]ΦTδ (k − 1, i + 1)ATδ [k − 1]l 1/2
− l, A[k − 1]Φδ (k − 1, i + 1)Rα [i]ΦTδ (k − 1, i + 1)AT [k − 1]l 1/2 )
+ l, Aδ [k − 1]Φδ (k − 1, k0 )X0 ΦTδ (k − 1, k0 )ATδ [k − 1]l 1/2
− l, A[k − 1]Φδ (k − 1, k0 )X0 ΦTδ (k − 1, k0 )AT [k − 1]l 1/2 .
Matrices Φδ (k −1, i+1)Rα [i]ΦTδ (k −1, i+1), k0 ≤ i < k −1, and Φδ (k −1, k0 )X0 ΦTδ (k −1, k0 ) are symmetric positive definite. We use lemma 6.1: for any ε1 > 0 and k, k ≥ k0 , there exists δk−1 such that −ε1 < ρ(l | Aδ [k − 1]Xδ,α (k − 1, k0 , E(x0 , X0 ))) − ρ(l | A[k − 1]Xδ (k − 1, k0 , E(x0 , X0 ))) < ε1 . If we denote the maximum of the largest singular values of matrices Φδ (k − 1, i + 1)Rα [i]ΦTδ (k − 1, i + 1), k0 ≤ i < k − 1, and
Φδ (k − 1, k0 )X0 ΦTδ (k − 1, k0 )
by hk−1 , and the maximum of the largest singular values of the matrices T T Uk−1 +Uk−1 Vk−1 Φδ (k−1, i+1)Rα [i]ΦTδ (k−1, i+1)AT [k−1], A[k−1]Φδ (k−1, i+1)Rα [i]ΦTδ (k−1, i+1)Vk−1
k0 ≤ i < k − 1, and T T Uk−1 +Uk−1 Vk−1 Φδ (k−1, k0 )X0 ΦTδ (k−1, k0 )AT [k−1] A[k−1]Φδ (k−1, k0 )X0 ΦTδ (k−1, k0 )Vk−1
by λk−1 , we have δk−1 =
−λk−1 +
2
ε1 2 λk−1 + 4( k−k ) 0
2hk−1 + 1
.
(38)
Remark. If we could find upper bounds h ≥ hi ,
λ ≥ λi ,
for all i, k0 ≤ i < k, δ would not depend on i: δk0 = · · · = δk−1 = δ =
−λ +
2
ε1 2 λ + 4( k−k ) 0
2h + 1
.
(39)
We now consider the third component of the difference (37). Let µ(R[i]) be the largest singular value of the matrix R[i], and define µmax = max {µ(R[i])}. k0 ≤i 0 there exist sequence {δi > 0}, k0 ≤ i < k, and parameter α > 0 such that 0 ≤ ρ(l | Xδ,α (k, k0 , E(x0 , X0 ))) − ρ(l | X (k, k0 , E(x0 , X0 ))) < ε for all vectors l and k ≥ k0 . Values for δi can be found from (38), α is defined in (41), with ε1 determined from the inequality
ε1 +
ε21 + 2µmax ε1
0. The corresponding control and the initial condition, which drive the system along the good trajectory are given by (25) and (26), with Φ substituted by Φδ and R by Rα . Formula (26) remains unchanged, and (25) is modified as uδ,α [i] =
7
P [i]B T [i]ΦTδ (k0 , i + 1)l0 + p[i]. l0 , Φδ (k0 , i + 1)Rα [i]ΦTδ (k0 , i + 1)l0 1/2
(47)
Example: regularized reach set
Consider the system
x1 [k + 1] x2 [k + 1]
=
0 1 0 0
x1 [k] x2 [k]
+
0 1
k ≥ 0,
u[k],
(48)
where |u[k]| ≤ µ for all k, and x[0] ∈ E(0, I). The matrices R[i] are now constant and have the form 0 0 . R[i] = R = 0 µ Since Ak = 0, k > 1, the reach set X (k, 0, E(0, I)) is a sum of two degenerate ellipsoids, X (k, 0, E(0, I)) = E(0, ARAT ) + E(0, R),
with T
ARA =
0 1 0 0
0 0 0 µ
0 0 1 0
=
µ 0 0 0
.
As seen in figure 2 (a), the reach set X (k, 0, E(0, I)) is actually a square with side 2µ and centered at the origin. The singular value decomposition of A is
A = U ΣV =
1 0 0 1
1 0 0 0
0 1 −1 0
Thus, by (33), Aδ can be expressed as
Aδ = A + δU V =
16
0 1+δ −δ 0
.
.
(a)
(b)
2.5 2 1.5 1 x 2
0.5 0
0. 5 1 1. 5 2 2. 5 10
5
0 x 1
5
10
Figure 2: (a) Reach set X (k, 0, E(0, I)) of the system (48) for k > 1 (µ = 1). (b) Regularized reach set Xδ,α at k = 10 with ε1 = 0.1. The expression for Rα comes from (27),
Rα = R + αI =
α 0 0 µ+α
.
The choice of α and δ depends on the desired accuracy ε and the number of steps k (in our example k0 = 0) we want this accuracy to be maintained. So, α is determined from (41), wherein µmax = µ. And δ is given by (39), wherein λ = 2µ and h = 1. In both formulas ε1 should satisfy condition (42). For example, for µ = 1, the choice of ε1 = 10−10 guarantees accuracy ε ≤ 0.01 of the exact reach set approximation after 100 steps. Figure 2 (b) shows the regularized reach set at time step k = 10. Here we picked ε1 = 0.1, calculated the corresponding δ = 0.000025 and α = 0.45 and used the recurrence relation (43) to calculate the shape matrix of the external ellipsoid for different direction vectors l0 .
8
Internal ellipsoidal approximations
Consider the sum of k + 1 ellipsoids E(0, Hi ), 0 ≤ i ≤ k with positive semidefinite Hi . An internal ellipsoidal approximation of this sum is given by E(0, Q) ⊆
k
E(0, Hi ),
i=0
with Q = Q[k] =
k i=0
k 1/2 T
Si H i
i=0
1/2
Si H i
and any orthogonal matrices Si , 0 ≤ i ≤ k, SiT Si = I (see [1]). Q[k] satisfies the recurrence relation Q[k + 1] =
k+1
k+1 1/2 T 1/2 Si H i Si H i i=0 i=0
17
,
(49)
k
=
i=0
1/2
Si H i
1/2
+ Sk+1 Hk+1
= Q[k] + Hk+1 +
+
1/2
Sk+1 Hk+1
T
k
i=0 k i=0
k T i=0
1/2 T
Si H i
1/2
1/2
Si H i
1/2
+ Sk+1 Hk+1
1/2
Sk+1 Hk+1
Si H i
,
(50)
with Q[0] = H0 . The following tightness condition is true (see [1]). Theorem 8.1 Let H0 be positive definite. Then the inequality
l, Q[k]l =
l,
k
k 1/2 T 1/2 Si H i Si H i l i=0 i=0
≤
k
l, Hi l 1/2
2
,
(51)
i=0
which holds for any l ∈ Rn and any array of orthogonal matrices Si , 0 ≤ i ≤ k, turns into equality for a given l ∈ Rn iff for all Hi there exist numbers ηi such that 1/2
Si H i
1/2
l = ηi S0 H0 l,
0 ≤ i ≤ k.
(52)
Now let us return to system (1) with nonsingular matrices A[i], k0 ≤ i < k, and its reach set at time step k, X (k, k0 , E(x0 , X0 )) given by (6). As we can see, the reach set is a sum of k + 1 ellipsoids. If we define Q[k] as
1/2
Q[k] = Φ(k, k0 ) X0 S0T + ×
1/2
S0 X0
+
k−1
k−1 i=k0
(Φ(k0 , i + 1)R[i]ΦT (k0 , i + 1))1/2 SiT [k]
Si [k](Φ(k0 , i + 1)R[i]ΦT (k0 , i + 1))1/2 ΦT (k, k0 ),
(53)
i=k0
with Q[k0 ] = X0 , and in which S0 and Si [k], k0 ≤ i < k, are any orthogonal matrices, E(q[k], Q[k]) ⊆ E(Φ(k, k0 )x0 , Φ(k, k0 )X0 ΦT (k, k0 )) +E
k−1
Φ(k, i + 1)B[i]p[i],
i=k0
Φ(k, i + 1)R[i]ΦT (k, i + 1)
i=k0
= X (k, k0 , E(x0 , X0 )), with q[k] coming from (7). X (k, k0 , E(x0 , X0 )).
k−1
(54)
So, the ellipsoid E(q[k], Q[k]) lies inside the reach set
18
Definition 8.1 Ellipsoid E(q[k], Q[k]) ⊆ X (k, k0 , E(x0 , X0 )) is a tight internal approximation of the reach set X (k, k0 , E(x0 , X0 )) if there exists direction l such that ρ(±l | E(q[k], Q[k])) = ρ(±l | X (k, k0 , E(x0 , X0 ))). In order for E(q[k], Q[k]) to be a tight internal approximation of the reach set for a given direction l, we use Theorem 8.1 and come up with the following tightness condition that can be checked by direct calculation. Theorem 8.2 For a given time step k ≥ k0 and given direction l ρ(l | E(q[k], Q[k])) = ρ(l | X (k, k0 , E(x0 , X0 ))) iff S0 and Si [k], k0 ≤ i < k, satisfy the relation 1/2
Si [k](Φ(k0 , i + 1)R[i]ΦT (k0 , i + 1))1/2 ΦT (k, k0 )l = ηi [k]S0 X0 ΦT (k, k0 )l, with ηi [k] =
l, Φ(k, i + 1)R[i]ΦT (k, i + 1)l 1/2 . l, Φ(k, k0 )X0 ΦT (k, k0 )l 1/2
(55)
(56)
Let l[k] satisfy the good curve relation (17). In this case, expressions (55) and (56) simplify into 1/2 (57) Si (Φ(k0 , i + 1)R[i]ΦT (k0 , i + 1))1/2 l0 = ηi S0 X0 l0 , and ηi =
l0 , Φ(k0 , i + 1)R[i]ΦT (k0 , i + 1)l0 1/2 . l0 , X0 l0 1/2
(58)
Now ηi and Si , k0 ≤ i < k no longer depend on k. We are ready to write the recurrence relation for the shape matrix Q[k] of the internal ellipsoid E(q[k], Q[k] that touches the reach set X (k, k0 , E(x0 , X0 )) at points of the good trajectory specified by (17). Using definition (53) and relation (50) we get Q[k + 1] = A[k]Q[k]AT [k] + R[k]
+ Φ(k + 1, k0 ) Mk Sk (Φ(k0 , k + 1)R[k]ΦT (k0 , k + 1))1/2
+ (Φ(k0 , k + 1)R[k]ΦT (k0 , k + 1))1/2 SkT MkT ΦT (k + 1, k0 ),
(59)
with Q[k0 ] = X0 , and 1/2
Mk = X0 S0T +
k−1 i=k0
(Φ(k0 , i + 1)R[i]ΦT (k0 , i + 1))1/2 SiT .
(60)
The orthogonal matrices S0 , Sk , k ≥ k0 , are related by (57) and (58). The recurrence relation for the center of the internal approximating ellipsoid q[k] is the same as for external and is given by (22). 19
Theorem 8.3 Let l[k] satisfy the equation (17), then the ellipsoid E(q[k], Q[k]), where Q[k] and q[k] are the solutions of (59) and (22), is such that for all k ≥ k0 , the following is true: E(q[k], Q[k]) ⊆ X (k, k0 , E(x0 , X0 )) and ρ(l[k] | X (k, k0 , E(x0 , X0 ))) − ρ(l[k] | E(q[k], Q[k])) = 0. The control that drives the system along the good trajectory is indicated in (25) for those vectors l0 , for which l0 , Φ(k0 , i + 1)R[i]ΦT (k0 , i + 1)l0 = 0,
k0 ≤ i < k,
and the initial condition in (26). For every boundary point of the reach set X (k, k0 , E(x0 , X0 )) there exists direction l[k] and corresponding tight internal ellipsoid E(q[k], Q[k]). Thus, we arrive at the result similar to theorem 4.2. Theorem 8.4 The reach set X (k, k0 , E(x0 , X0 )) can be expressed as X (k, k0 , E(x0 , X0 )) =
E(q[k], Q[k]).
(61)
l[k]
Remark. The results of this section are also valid for singular matrices R[i], and no special regularization as in the case of external approximation is needed. Until now, in discussing internal ellipsoidal approximations, we assumed A[i] to be nonsingular. In case of singular A[i], the recurrence relations above lose their sense because matrix Φ(k0 , k) does not exist. It is still possible to build internal approximating ellipsoids using (49) even for the singular case, since the reach set of a linear discrete-time system for every time step k is nothing else but a finite sum of ellipsoids. However, the recurrence relations in the case of singular matrix A cease to exist, making the calculation of tight internal approximating ellipsoids a rather involved procedure because for every k the matrices Si need to be recomputed. Therefore, we suggest the following way out for singular case. Instead of the actual reach set X (k, k0 , E(x0 , X0 )) we consider the regularized reach set Xδ,α (k, k0 , E(x0 , X0 )) defined in (35) and constructed according to theorem 6.1. Matrices Aδ [i] are nonsingular, so are Rα [i], and the relations (57)-(60). Recall the example (48) from section 7. Figure 3 illustrates the internal ellipsoidal approximations √ ofthe regularized 1 3 1 2 √ 2 , and . reach set at k = 10 for three different initial directions, l0 = 3 0 − 12 2
20
1 0.8 0.6 0.4
x2
0.2 0 −0.2 −0.4 −0.6 −0.8 −1
−4
−3
−2
−1
0 x
1
2
3
4
1
Figure 3: Internal ellipsoidal approximations of the regularized reach set Xδ,α at k = 10.
9
Conclusion
This paper extends the results of [1] to linear discrete-time systems. We do not assume controllability of the system and also consider the situation when the state transition matrix Φ(k, k0 ) is singular. It was shown that in these situations it is possible to approximate the actual reach set with any given ε-accuracy by regularizing matrices B[i]P [i]B T [i] and A[i]. Theorem 6.1 provides the means of picking the appropriate regularization parameters. Internal ellipsoidal approximations work well even with singular matrices B[i]P [i]B T [i] and Φ(k, k0 ). However, in order to use the concept of ‘good curves’ and the recurrence relations (57)-(60), we still need to perform the regularization.
References [1] A. Kurzhanski and P. Varaiya, “On ellipsoidal techniques for reachability analysis,” Optimization Methods and Software, vol. 17, pp. 177–237, 2000. [2] N. Giorgetti and G. Pappas, “Bounded model checking for hybrid dynamical systems,” in IEEE Conference on Decision and Control, Seville, Spain, 2005, submitted. [3] I. Mitchell and C. Tomlin, “Level set methods for computation in hybrid systems,” in Hybrid Systems: Computation and Control, ser. Lecture Notes in Computer Science, N. Lynch and B. Krogh, Eds. Springer, 2000, vol. 1790, pp. 21–31. [4] “Level set toolbox,” www.cs.ubc.ca/˜mitchell/ToolboxLS.
21
[5] M. Kvasnica, P. Grieder, M. Baoti´c, and M. Morari, “Multi parametric toolbox (mpt),” in Hybrid Systems: Computation and Control, R. Alur and G. Pappas, Eds. Springer, 2004, vol. 2993, pp. 448–462. [6] “Mpt homepage,” control.ee.ethz.ch/˜mpt. [7] T. Motzkin, H. Raiffa, G. Thompson, and R. Thrall, “The double description method,” in Conttributions to Theory of Games, H.W.Kuhn and A.W.Tucker, Eds. Princeton University Press, 1953, vol. 2. [8] “Cdd/cdd+ homepage,” www.cs.mcgill.ca/˜fukuda/soft/cdd home/cdd.html. [9] “d/dt homepage,” www-verimag.imag.fr/˜tdang/ddt.html. [10] E. Asarin, O. Bournez, T. Dang, and O. Maler, “Approximate reachability analysis of piecewise linear dynamical systems,” in Hybrid Systems: Computation and Control, ser. Lecture Notes in Computer Science, N.Lynch and B.H.Krogh, Eds. Springer, 2000, vol. 1790, pp. 21–31. [11] D. Avis, D. Bremner, and R. Seidel, “How good are convex hull algorithms?” Computational Geometry: Theory and Applications, vol. 7, pp. 265–301, 1997. [12] A. Girard, “Reachability of uncertain linear systems using zonotopes,” in Hybrid Systems: Computation and Control, ser. Lecture Notes in Computer Science. Springer, 2005, to appear. [13] “Zonotope methods on wolfgang k¨ uhn homepage,” www.decatur.de. [14] “CheckMate homepage,” www.ece.cmu.edu/˜webk/checkmate. [15] O. Stursberg and B. Krogh, “Efficient representation and computation of reachable sets for hybrid systems,” in Hybrid Systems: Computation and Control, ser. Lecture Notes in Computer Science, O.Maler and A.Pnueli, Eds. Springer, 2003, vol. 2623, pp. 482–497. [16] “Requiem homepage,” www.seas.upenn.edu/˜hybrid/requiem/requiem.html. [17] G. Lafferriere, G. Pappas, and S. Yovine, “Symbolic reachability computation for families of linear vector fields,” Journal of Symbolic Computation, vol. 32, pp. 231–253, 2001. [18] E. Kostousova, “Control synthesis via parallelotopes: optimization and parallel computations,” Optimization Methods and Software, vol. 14, no. 4, pp. 267–310, 2001. [19] P. Varaiya, “Reach set computation using optimal control,” Proc. of KIT Workshop on Verification on Hybrid Systems. Verimag, Grenoble., 1998. [20] A. Kurzhanski and P. Varaiya, “Dynamic optimization for reachability related problems,” Journal of Optimization, Theory and Applications, vol. 108, no. 2, pp. 227–251, 2001.
22
Abstract This paper describes the computation of reach sets for discrete-time linear control systems with time-varying coefficients and ellipsoidal bounds on the controls and initial conditions. The algorithms construct external and internal ellipsoidal approximations so that they touch the reach set boundary from outside and from inside. Recurrence relations that describe the time evolution of these approximations are provided. The paper also deals with discrete-time linear systems with singular state transition matrix.
1
Introduction
We extend the use of ellipsoidal methods for continuous-time linear systems developed in [1] to calculate the reach set for discrete-time linear systems. We also deal with systems with a singular state transition matrix and show how to approximate the reach set with any given accuracy using regularization and provide the means for picking the appropriate regularization parameters that guarantee this accuracy over a specified finite time interval. We make no controllability assumptions regarding the system. Thus, the suggested regularization technique also extends results of [1] to systems that are not completely controllable. In case of singular state transition matrix, or if the system is not controllable at every time step, the exact representation of the reach set on every time step by external ellipsoidal approximation is not possible. But the exact representation is possible for the regularized reach set that overapproximates the actual one, and whose boundary can be made arbitrarily close to the boundary of the actual one over any given finite time interval. Using internal ellipsoidal approximations, on the other hand, it is always possible to exactly represent ∗
Research supported by NSF Grant CCR-00225610.
1
the actual reach set. However, the concept of ‘good curves’ that was introduced in [1] and without which the computation process becomes heavy and complicated, can only be applied for systems with nonsingular state transition matrix. Therefore, in order to make use of good curves in the singular case, we perform the internal approximation of the regularized reach set.
2
Related work
Research on dynamical and hybrid systems has produced methods for verification and controller synthesis. A common step in most methods is the calculation of reach sets. For example, [2] uses bounded model checking for verification of discrete-time linear hybrid systems. For that, the hybrid system is converted into a discrete transition system. The conversion requires the computation of reach sets of discrete-time linear systems. Exact computation of reach sets is not generally possible, and approximation techniques are needed. The choice of the reach set representation determines the efficiency of such techniques. The more complex the representation is, the more costly is the storage of the sets, the more elaborate is the computation, but the better is the reach set approximation. Choosing the reach set representation is a compromise between these factors. The level set method [3, 4] deals with general nonlinear controlled systems and gives exact representation of their reach sets, but the computation is complex, requiring solving the HJB equation and finding the set of states that belong to sub-zero level set of the value function. The method [4] is impractical for systems of dimension higher than three. Reachability analysis of controlled linear systems with convex bounds on the control and initial conditions reduces to effective manipulation with convex sets, performing set-valued operations such as unions, intersections, geometric sums and differences. Two basic objects are used as convex approximations: various kinds of polytopes, e.g. general polytopes, zonotopes, parallelotopes, rectangular polytopes; and ellipsoids. Reachability analysis and verification for discrete-time linear systems through general polytopes are implemented in the Multi Parametric Toolbox (MPT) for Matlab [5], [6]. The reach set at every time step is computed as the geometric sum of two polytopes. The procedure consists in finding the vertices of the resulting polytope and calculating its convex hull. The convex hull algorithm employed by MPT is based on the Double Description method [7] and implemented in the CDD/CDD+ package [8]. Its complexity is V n , where V is the number of vertices and n is the state space dimension. Hence the use of MPT is practical for low dimensional systems. But even in low dimensional systems the number of vertices in the the reach set polytope can grow very large with the number of time steps. For example, consider the system, xk+1 = Axk + uk ,
2
cos 1 − sin 1 with A = , uk ∈ {u ∈ R2 | u∞ ≤ 1}, and x0 ∈ {x ∈ R2 | x∞ ≤ 1}. sin 1 cos 1 Starting with a rectangular initial set, the number of vertices of the reach set polytope is 4k + 4 at the kth step. In d/dt [9], the reach set is approximated by unions of rectangular polytopes [10]. The x 1(t 2)
R (t 2) x 2(t 2) R [t 1,t 2] R (t 1) x 1(t 1)
x 2(t 1)
x 1(t 1)
x 2(t 1)
R [t 0,t 1] x 1(t 0)
x 2(t 0)
(a)
x 1(t 0)
x 2(t 0)
(b)
(c)
(d)
(e)
(f)
Figure 1: Reach set approximation by union of rectangles. algorithm works as follows. First, given the set of initial conditions defined as a polytope, the evolution in time of the polytope’s extreme points is computed (figure 1(a)). R(t1 ) in figure 1(a) is the reach set of the system at time t1 , and R[t0 , t1 ] is the set of all points that can be reached during [t0 , t1 ]. Second, the algorithm computes the convex hull of vertices of both, the initial polytope and R(t1 ) (figure 1(b)). The resulting polytope is then bloated to include all the reachable states in [t0 , t1 ] (figure 1(c)). Finally, this overapproximating polytope is in its turn overapproximated by the union of rectangles (figure 1(d)). The same procedure is repeated for the next time interval [t1 , t2 ], and the union of both rectangular approximations is taken (figure 1(e,f)), and so on. Rectangular polytopes are easy to represent and the number of facets grows linearly with dimension, but a large number of rectangles must be used to assure the approximation is not overly conservative. Besides, the important part of this method is again the convex hull calculation whose implementation relies on the same CDD/CDD+ library. This limits the dimension of the system and time interval for which it is practical to calculate the reach set. Polytopes can give arbitrarily close approximations to convex sets, but the number of vertices can grow prohibitively large and, as shown in [11], the computation of a polytope by its convex hull becomes intractable for large number of vertices in high dimensions. The method of zonotopes for external approximation of reach sets is presented in [12]. A 3
zonotope is a special class of polytopes (see [13]) of the form, n
Z = {x ∈ R | x = c +
p
αi gi , − 1 ≤ αi ≤ 1},
i=1
wherein c and g1 , ..., gp are vectors in Rn . Thus, a zonotope Z is defined by its center c and generator vectors g1 , ..., gp . The value p/n is called the order of the zonotope. This compact representation and the fact that they are closed under linear transformation and geometric sum, make zonotopes attractive for reach set computation. The difficulty is that with every integration step the order of the approximating zonotope increases by p/n. The difficulty can be averted by limiting the number of generator vectors, and overapproximating zonotopes whose number of generator vectors exceeds this limit by lower order zonotopes. The zonotope method of [12] provides neither a notion of tightness of approximation nor an estimate of the error resulting from the zonotope order reduction. Further, the method does not indicate any way of control synthesis, i.e. finding a sequence of controls that drives the system to the desired terminal set or point. CheckMate [14] is a Matlab toolbox that can evaluate specifications for trajectories starting from the set of initial (continuous) states corresponding to the parameter values at the vertices of the parameter set. This provides preliminary insight into whether the specifications will be true for all parameter values. The method of oriented rectangluar polytopes for external approximation of reach sets is introduced in [15]. The basic idea is to construct an oriented rectangular hull of the reach set for every time step, whose orientation is determined by the singular value decomposition of the sample covariance matrix for the states reachable from the vertices of the initial polytope. The limitation of CheckMate and the method of oriented rectangles is that only autonomous (i.e. uncontrolled) systems are allowed, and only an external approximation of the reach set is provided. Requiem [16] is a Mathematica notebook which, given a linear system, the set of initial conditions and control bounds, symbolically computes the reach set exactly, using the experimental quantifier elimination package. Quantifier elimination is the removal of all quantifiers (the universal quantifier ∀ and the existential quantifier ∃) from a quantified system. Each quantified formula is substituted with quantifier-free expression with operations +, ×, = and k0 from the initial position (k0 , x0 ) is the set of all states x[k] reachable at time k by the system (1) with x[k0 ] = x0 through all possible controls that satisfy (3). The reach set X (k, k0 , X0 ) is X (k, k0 , X0 ) =
{X (k, k0 , x0 ) | x0 ∈ X0 )}.
Some important facts about reach sets follow. Lemma 3.1 The set-valued map X (k, k0 , X0 ) satisfies the semigroup property X (k, k0 , X0 ) = X (k, i, X (i, k0 , X0 )),
k0 ≤ i ≤ k.
(5)
Lemma 3.2 The reach set X (k, k0 , E(x0 , X0 )) can be expressed as X (k, k0 , E(x0 , X0 )) = q[k] + Φ(k, k0 )E(0, X0 ) +
k−1
Φ(k, i + 1)E(0, B[i]P [i]B T [i]),
(6)
i=k0
where q[k] = Φ(k, k0 )x0 +
k−1
Φ(k, i + 1)B[i]p[i].
(7)
i=k0
Using (6) and (4), we obtain the next result. Lemma 3.3 The support function of the reach set is ρ(l | X (k, k0 , E(x0 , X0 ))) = l, q[k] + l, Φ(k, k0 )X0 ΦT (k, k0 )l 1/2 +
k−1
l, Φ(k, i + 1)B[i]P [i]B T [i]ΦT (k, i + 1)l 1/2 .
i=k0
Representation (6) leads to the next fact. 6
(8)
Lemma 3.4 The reach set X (k, k0 , E(x0 , X0 )) is a convex compact set in Rn . If the matrices A[i], k0 ≤ i < k, are nonsingular, the ellipsoid Φ(k, k0 )E(0, X0 ) in sum (6) is nondegenerate, and so the reach set X (k, k0 , E(x0 , X0 )) contains an open set of Rn . Its boundary points have an important characterization: y is a boundary point of X (k, k0 , E(x0 , X0 )) only if there exists a support vector ly = 0 so that ly , y = ρ(ly | X (k, k0 , E(x0 , X0 ))) = max{ ly , x | x ∈ X (k, k0 , E(x0 , X0 ))}.
(9)
The control uy [i], k0 ≤ i < k, and the initial state x[k0 ] = x0y ∈ E(x0 , X0 ), which drive the system (1) from the state x[k0 ] = x0y to a boundary point x[k] = y is determined from the following theorem. Theorem 3.1 Suppose A[i], k0 ≤ i < k, are nonsingular, and x[k] = y is a boundary point of X (k, k0 , E(x0 , X0 )). Then the control uy [i], k0 ≤ i < k, and the initial state x[k0 ] = x0y , which yield the unique trajectory xy [i] that reaches y = xy [k] from x0y = x[k0 ], and ensures (9), satisfy the following maximum principle for control: ly , Φ(k, i + 1)B[i]uy [i] = max{ ly , Φ(k, i + 1)B[i]u | u ∈ Φ(k, i + 1)E(p[i], P [i])} = ly , Φ(k, i + 1)B[i]p[i]
+ ly , Φ(k, i + 1)B[i]P [i]B T [i]ΦT (k, i + 1)ly 1/2
(10)
for k0 ≤ i < k, and the maximum condition for the initial state: ly , Φ(k, k0 )x0y = max{ ly , x | x ∈ Φ(k, k0 )E(x0 , X0 )} = ly , Φ(k, k0 )x0 + ly , Φ(k, k0 )X0 ΦT (k, k0 )ly 1/2 .
(11)
Here ly is the support vector for the set X (k, k0 , E(x0 , X0 )) at the point y that satisfies (9). Remark. To simplify the notation we shall denote R[i] = B[i]P [i]B T [i] and from now on use R[i] instead of B[i]P [i]B T [i] where possible.
4
External ellipsoidal approximations
Relation (6) shows that the reach set X (k, k0 , E(x0 , X0 )) is a sum of k − k0 + 1 ellipsoids. We overapproximate the reach set by a single ellipsoid. Here we consider the situation when matrices A[i] are nonsingular. The case with singular A[i] is treated in section 6. In the beginning we shall also assume that matrices R[i] are nonsingular. Remark. R[i] are nonsingular if the matrices B[i] ∈ Rn×m have rank n (n ≤ m). We shall look for external approximating ellipsoids for the reach set X (k, k0 , E(x0 , X0 )) in the class of ellipsoids E(q[k], Q[k]) that satisfy the following two relations: q[k] = A[k − 1]q[k − 1] + B[k − 1]p[k − 1], 7
q[k0 ] = x0 ,
(12)
and Q[k] =
k−1
si [k] + s0 [k]
i=k0
×
k−1 i=k0
T s−1 i [k]Φ(k, i + 1)R[i]Φ (k, i + 1)
T +s−1 0 [k]Φ(k, k0 )X0 Φ (k, k0 ) ,
Q[k0 ] = X0 ,
(13)
for arbitrary s0 [k] > 0, si [k] > 0 for k > k0 and k0 ≤ i < k. Using (8), it can be easily checked that ρ(l | E(q[k], Q[k])) ≥ ρ(l | X (k, k0 , E(x0 , X0 ))) for all directions l. Hence, X (k, k0 , E(x0 , X0 )) ⊆ E(q[k], Q[k]). Definition 4.1 We say that ellipsoid E(q[k], Q[k]), defined by (12) and (13), is a tight external approximation of the reach set X (k, k0 , E(x0 , X0 )) if there exists direction l such that ρ(±l | E(q[k], Q[k])) = ρ(±l | X (k, k0 , E(x0 , X0 ))). A tight external ellipsoid defined by direction l, touches the boundary of the reach set at point xl such that ρ(l | X (k, k0 , E(x0 , X0 ))) = l, xl = ρ(l | E(q[k], Q[k])) = l, q[k] + l, Q[k]l 1/2 . This condition will be fulfilled if we choose
and
si [k] = l, Φ(k, i + 1)R[i]ΦT (k, i + 1)l 1/2 ,
(14)
s0 [k] = l, Φ(k, k0 )X0 ΦT (k, k0 )l 1/2 .
(15)
To check it, we substitute (14) and (15) into (13), calculate the support function of the resulting ellipsoid and arrive at (8). As we can see, for every direction l there exists a tight external ellipsoidal approximation of the reach set X (k, k0 , E(x0 , X0 )), whose shape matrix is calculated by (13). For all 0 ≤ i < k, the parameters si [k] in (15) and (14) depend on k. As a result, the sum (13) needs to be recomputed for every k with new values of s0 [k] and si [k]. To get rid of this dependency on k, and thus greatly reduce the computation effort, we introduce the notion of ‘good curves’. These are trajectories l[k] that satisfy the homogeneous adjoint equation l[k] = AT [k]l[k + 1], 8
l(k0 ) = l0 .
(16)
Vector l[k] can be expressed as l[k] = ΦT (k0 , k)l0
(for time-invariant systems: l[k] = ((Ak−k0 )T )−1 l0 ).
(17)
With l[k] defined in (17), relations (14) and (15) simplify into si = l0 , Φ(k0 , i + 1)R[i]Φ(k0 , i + 1)l0 1/2
(18)
s0 = l0 , X0 l0 1/2 .
(19)
and
Now si and s0 no longer depend on k. When substituted into (13), they give us a recurrence relation for matrix Q[k], Q[k + 1] = (1 + π[k])A[k]Q[k]AT [k] + (1 + π −1 [k])R[k],
Q[k0 ] = X0 ,
(20)
in which π[k] is defined by sk i=k0 si + s0
π[k] =
k−1
=
k−1
=
l0 , Φ(k0 , k + 1)R[k]ΦT (k0 , k + 1)l0 1/2 . l0 , Φ(k0 , k)Q[k]ΦT (k0 , k)l0 1/2
l0 , Φ(k0 , k + 1)R[k]ΦT (k0 , k + 1)l0 1/2 T 1/2 + l , X l 1/2 0 0 0 i=k0 l0 , Φ(k0 , i + 1)R[i]Φ (k0 , i + 1)l0
(21)
The expression for the center of the approximating ellipsoid is given by q[k + 1] = A[k]q[k] + B[k]p[k],
q[k0 ] = x0 .
(22)
In both (20) and (22) k ≥ k0 . We state this result as a theorem. Theorem 4.1 Let l[k] satisfy the equation (17). Then the ellipsoid E(q[k], Q[k]) given by (20) and (22), is such that for all k ≥ k0 the reach set X (k, k0 , E(x0 , X0 )) ⊆ E(q[k], Q[k]) and ρ(l[k] | E(q[k], Q[k])) − ρ(l[k] | X (k, k0 , E(x0 , X0 ))) = 0. The supporting hyperplane for the reach set X (k, k0 , E(x0 , X0 )) generated by the vector l[k] is also a supporting hyperplane for E(q[k], Q[k]) and touches it at the same point xl [k], given by Q[k]l[k] . (23) xl [k] = q[k] + l[k], Q[k]l[k] 1/2 Since for every boundary point xl [k] of the reach set X (k, k0 , E(x0 , X0 )) there exists direction l[k] and corresponding tight external ellipsoid E(q[k], Q[k]), the following result is established. 9
Theorem 4.2 The reach set X (k, k0 , E(x0 , X0 )) can be expressed as X (k, k0 , E(x0 , X0 )) =
E(q[k], Q[k]).
(24)
l[k]
If l[k] is a “good curve”, it specifies the external ellipsoid E(q[k], Q[k]), with Q[k] and q[k] coming from (20) and (22), which touches the reach set at some point xl [k]. In order for the points xl [·] to form the system trajectory on k0 , ..., k, that is, for xl [k] to satisfy xl [k] = Φ(k, k0 )x0 +
k−1
Φ(k, i + 1)B[i]u[i],
i=k0
the corresponding control u[i] ∈ E(p[i], P [i]) and the initial condition x0 ∈ E(x0 , X0 ) are derived from the maximum principle for control (10) and maximum condition (11). They are specified by P [i]B T [i]ΦT (k0 , i + 1)l0 + p[i], (25) u[i] = l0 , Φ(k0 , i + 1)R[i]ΦT (k0 , i + 1)l0 1/2 and x0l =
X0 l0 + x0 . l0 , X0 l0 1/2
(26)
Thus, if the system starts at x0l at time step k0 , the control u[i] will bring it at time step k to the boundary point of the reach set, at which the reach set is touched by the external approximating ellipsoid E(q[k], Q[k]) defined by the direction l[k]. Remark. For n-dimensional system, every step of computation of external ellipsoidal approximation requires 2n2 scalar multiplications for the center of ellipsoid and (8n3 + 4n2 + 2n) scalar multiplications for its shape matrix. Since the trajectory of the center is the same for different values of parameter l0 , computing L ellipsoidal approximations for k time steps requires k(L(8n3 + 4n2 + 2n) + 2n2 ) scalar multiplications. So, the complexity of computation grows linearly with the number of time steps and number of approximations, and polynomially, as O(n3 ), with the state space dimension.
5
Singular R
We start by stating a simple fact that will be useful in our proofs. Lemma 5.1 Let a ≥ b > 0. Then from a2 − b2 < ε2 for some ε > 0, it follows that a − b < ε. 10
In case the matrix R[i] is singular for some i, the recurrence relations (20) are not valid because we cannot guarantee that parameter si defined in (18) is strictly positive. Thus, we need to replace R[i] with a suitable nonsingular matrix. Define Rα [i], Rα [i] = R[i] + α2 I,
k0 ≤ i < k,
(27)
for α > 0. Rα [i] is positive definite and Rα [i] → R[i] as α → 0. Corresponding to Rα [i] we define the “regularized” reach set Xα (k, k0 , E(x0 , X0 )), Xα (k, k0 , E(x0 , X0 )) = q[k] + Φ(k, k0 )E(0, X0 ) +
k−1
Φ(k, i + 1)E(0, Rα [i]),
(28)
i=k0
with support function ρ(l | Xα (k, k0 , E(x0 , X0 ))) = l, q[k] + l, Φ(k, k0 )X0 ΦT (k, k0 )l 1/2 +
k−1
l, Φ(k, i + 1)Rα [i]ΦT (k, i + 1)l 1/2 ,
(29)
i=k0
with q[k] determined from (7). The difference between the support functions of the regularized reach set Xα (k, k0 , E(x0 , X0 )) and the actual reach set X (k, k0 , E(x0 , X0 )) is ρ(l | Xα (k, k0 , E(x0 , X0 ))) − ρ(l | X (k, k0 , E(x0 , X0 ))) =
k−1
( l, Φ(k, i + 1)Rα [i]ΦT (k, i + 1)l 1/2 − l, Φ(k, i + 1)R[i]ΦT (k, i + 1)l 1/2 ).
i=k0
Let us estimate the ith element of this sum, l, Φ(k, i + 1)Rα [i]ΦT (k, i + 1)l 1/2 − l, Φ(k, i + 1)R[i]ΦT (k, i + 1)l 1/2 ,
k0 ≤ i < k.
Let σ(A[i]) be the largest singular value of matrix A[i], and denote σmax = max {σ(A[i])}. k0 ≤i 0 there exists α > 0, such that 0 ≤ ρ(l | Xα (k, k0 , E(x0 , X0 ))) − ρ(l | X (k, k0 , E(x0 , X0 ))) < ε for all k ≥ k0 , and all vectors l. More precisely, α can be chosen as ε , α= k−k0 −1 2 2(σmax + · · · + σmax + σmax + 1) where σmax is defined in (30). 11
(31)
Remark. In (30) we could have taken i strictly greater than k0 . We leave this definition as is, however, because it will be used in the next section in its present form. Substituting R with Rα in (18) and (21), we obtain recurrence relations for Qα [k], the shape matrix of the tight external ellipsoid of the regularized reach set Xα (k, k0 , E(x0 , X0 )). The recurrence relation for the center of this ellipsoid (22) remains untouched. So, the relation X (k, k0 , E(x0 , X0 )) =
E(q[k], Qα [k])
(32)
l[k]
holds with ε accuracy in the Hausdorff metric. Expression (25) for the control is similarly modified by replacing R with Rα .
6
Singular A
When A[i] is singular for some i, k0 ≤ i < k, relation (17) is invalid because the inverse for the state transition matrix Φ(k, k0 ) does not exist; hence, we cannot use the expression Φ(k0 , k); and the subsequent formulas (18) and (20) do not make sense for singular A[i]. The way around it is to regularize matrix A[i]. Let the singular value decomposition of A[i] be A[i] = Ui Σi Vi , and define the new matrix Aδ [i] by Aδ [i] = Ui (Σi + δI)Vi = A[i] + δUi Vi ,
k0 ≤ i < k,
(33)
where δ > 0. As one can see, Aδ [i] is nonsingular for δ > 0, and Aδ [i] → A[i] as δ → 0. Corresponding to the matrix Aδ [i] there is an invertible state transition matrix Φδ (k, k0 ). Remark. Here, actually, δ may depend on i, and Φδ (k, k0 ) is Φδ (k, k0 ) = Aδk−1 [k − 1]Aδk−2 [k − 2] · · · Aδk0 +1 [k0 + 1]Aδk0 [k0 ]. We drop the subindices, however, to simplify the notation. Let H be any symmetric positive definite matrix. We would like to find such δ > 0 that for given ε > 0, −ε2 < l, Aδ [i]HATδ [i]l − l, A[i]HAT [i]l < ε2 . for all vectors l, and any i, then by lemma 5.1, −ε < l, Aδ [i]HATδ [i]l 1/2 − l, A[i]HAT [i]l 1/2 < ε, For that we write Aδ [i]HATδ [i] − A[i]HAT [i] = δ(A[i]HViT UiT + Ui Vi HAT [i]) + δ2 Ui Vi HViT UiT . 12
If we denote the largest singular value of H as h, and the largest singular value of the symmetric matrix A[i]HViT UiT + Ui Vi HAT [i] as λ = λ(A[i]HViT UiT + Ui Vi HAT [i]), then δ could be found from the inequality hδ2 + λδ < ε2 . Thus, we arrive at the following lemma. Lemma 6.1 Let H = H T > 0. Then, for any ε > 0 there exists δ > 0 such that −ε < l, Aδ [i]HATδ [i]l 1/2 − l, A[i]HAT [i]l 1/2 < ε. More precisely, δ can be chosen as
2
−λ + λ + 4ε2 , δ= 2h + 1
(34)
where λ is the largest singular value of the matrix A[i]HViT UiT + Ui Vi HAT [i]. We define the regularized reach set Xδ,α (k, k0 , E(x0 , X0 )) by Xδ,α(k, k0 , E(x0 , X0 )) = q[k] + Φδ (k, k0 )E(0, X0 ) +
k−1
Φδ (k, i + 1)E(0, Rα [i]).
(35)
i=k0
Its support function is ρ(l | Xδ,α (k, k0 , E(x0 , X0 ))) = l, q[k] + l, Φδ (k, k0 )X0 ΦTδ (k, k0 )l 1/2 +
k−1 i=k0
l, Φδ (k, i + 1)Rα [i]ΦTδ (k, i + 1)l 1/2 ,
(36)
where Rα [i] is defined in (27) and q[k] comes from (7) as before. Our goal is to make the regularized reach set arbitrary close to the actual reach set for every k ≥ k0 , by the appropriate choice of parameters δi , k0 ≤ i < k and α. In other words, for any ε > 0 we should find δi > 0, k0 ≤ i < k, and α > 0, such that 0 ≤ ρ(l | Xδ,α (k, k0 , E(x0 , X0 ))) − ρ(l | X (k, k0 , E(x0 , X0 ))) < ε for all k ≥ k0 , and all vectors l. We now show how to construct such α and the sequence δi . Let us estimate the difference of the support functions of the regularized reach set and the actual reach set ρ(l | Xδ,α (k, k0 , E(x0 , X0 ))) − ρ(l | X (k, k0 , E(x0 , X0 ))) = ρ(l | Aδ [k − 1]Xδ,α (k − 1, k0 , E(x0 , X0 ))) − ρ(l | A[k − 1]Xδ,α (k − 1, k0 , E(x0 , X0 ))) +ρ(l | A[k − 1]Xδ,α (k − 1, k0 , E(x0 , X0 ))) − ρ(l | A[k − 1]X (k − 1, k0 , E(x0 , X0 ))) +ρ(l | E(0, Rα [k − 1])) − ρ(l | E(0, R[k − 1])). 13
(37)
As we can see, this difference breaks up into three components. First, we analyze ρ(l | Aδ [k − 1]Xδ,α (k − 1, k0 , E(x0 , X0 ))) − ρ(l | A[k − 1]Xδ,α (k − 1, k0 , E(x0 , X0 ))) =
k−2 i=k0
( l, Aδ [k − 1]Φδ (k − 1, i + 1)Rα [i]ΦTδ (k − 1, i + 1)ATδ [k − 1]l 1/2
− l, A[k − 1]Φδ (k − 1, i + 1)Rα [i]ΦTδ (k − 1, i + 1)AT [k − 1]l 1/2 )
+ l, Aδ [k − 1]Φδ (k − 1, k0 )X0 ΦTδ (k − 1, k0 )ATδ [k − 1]l 1/2
− l, A[k − 1]Φδ (k − 1, k0 )X0 ΦTδ (k − 1, k0 )AT [k − 1]l 1/2 .
Matrices Φδ (k −1, i+1)Rα [i]ΦTδ (k −1, i+1), k0 ≤ i < k −1, and Φδ (k −1, k0 )X0 ΦTδ (k −1, k0 ) are symmetric positive definite. We use lemma 6.1: for any ε1 > 0 and k, k ≥ k0 , there exists δk−1 such that −ε1 < ρ(l | Aδ [k − 1]Xδ,α (k − 1, k0 , E(x0 , X0 ))) − ρ(l | A[k − 1]Xδ (k − 1, k0 , E(x0 , X0 ))) < ε1 . If we denote the maximum of the largest singular values of matrices Φδ (k − 1, i + 1)Rα [i]ΦTδ (k − 1, i + 1), k0 ≤ i < k − 1, and
Φδ (k − 1, k0 )X0 ΦTδ (k − 1, k0 )
by hk−1 , and the maximum of the largest singular values of the matrices T T Uk−1 +Uk−1 Vk−1 Φδ (k−1, i+1)Rα [i]ΦTδ (k−1, i+1)AT [k−1], A[k−1]Φδ (k−1, i+1)Rα [i]ΦTδ (k−1, i+1)Vk−1
k0 ≤ i < k − 1, and T T Uk−1 +Uk−1 Vk−1 Φδ (k−1, k0 )X0 ΦTδ (k−1, k0 )AT [k−1] A[k−1]Φδ (k−1, k0 )X0 ΦTδ (k−1, k0 )Vk−1
by λk−1 , we have δk−1 =
−λk−1 +
2
ε1 2 λk−1 + 4( k−k ) 0
2hk−1 + 1
.
(38)
Remark. If we could find upper bounds h ≥ hi ,
λ ≥ λi ,
for all i, k0 ≤ i < k, δ would not depend on i: δk0 = · · · = δk−1 = δ =
−λ +
2
ε1 2 λ + 4( k−k ) 0
2h + 1
.
(39)
We now consider the third component of the difference (37). Let µ(R[i]) be the largest singular value of the matrix R[i], and define µmax = max {µ(R[i])}. k0 ≤i 0 there exist sequence {δi > 0}, k0 ≤ i < k, and parameter α > 0 such that 0 ≤ ρ(l | Xδ,α (k, k0 , E(x0 , X0 ))) − ρ(l | X (k, k0 , E(x0 , X0 ))) < ε for all vectors l and k ≥ k0 . Values for δi can be found from (38), α is defined in (41), with ε1 determined from the inequality
ε1 +
ε21 + 2µmax ε1
0. The corresponding control and the initial condition, which drive the system along the good trajectory are given by (25) and (26), with Φ substituted by Φδ and R by Rα . Formula (26) remains unchanged, and (25) is modified as uδ,α [i] =
7
P [i]B T [i]ΦTδ (k0 , i + 1)l0 + p[i]. l0 , Φδ (k0 , i + 1)Rα [i]ΦTδ (k0 , i + 1)l0 1/2
(47)
Example: regularized reach set
Consider the system
x1 [k + 1] x2 [k + 1]
=
0 1 0 0
x1 [k] x2 [k]
+
0 1
k ≥ 0,
u[k],
(48)
where |u[k]| ≤ µ for all k, and x[0] ∈ E(0, I). The matrices R[i] are now constant and have the form 0 0 . R[i] = R = 0 µ Since Ak = 0, k > 1, the reach set X (k, 0, E(0, I)) is a sum of two degenerate ellipsoids, X (k, 0, E(0, I)) = E(0, ARAT ) + E(0, R),
with T
ARA =
0 1 0 0
0 0 0 µ
0 0 1 0
=
µ 0 0 0
.
As seen in figure 2 (a), the reach set X (k, 0, E(0, I)) is actually a square with side 2µ and centered at the origin. The singular value decomposition of A is
A = U ΣV =
1 0 0 1
1 0 0 0
0 1 −1 0
Thus, by (33), Aδ can be expressed as
Aδ = A + δU V =
16
0 1+δ −δ 0
.
.
(a)
(b)
2.5 2 1.5 1 x 2
0.5 0
0. 5 1 1. 5 2 2. 5 10
5
0 x 1
5
10
Figure 2: (a) Reach set X (k, 0, E(0, I)) of the system (48) for k > 1 (µ = 1). (b) Regularized reach set Xδ,α at k = 10 with ε1 = 0.1. The expression for Rα comes from (27),
Rα = R + αI =
α 0 0 µ+α
.
The choice of α and δ depends on the desired accuracy ε and the number of steps k (in our example k0 = 0) we want this accuracy to be maintained. So, α is determined from (41), wherein µmax = µ. And δ is given by (39), wherein λ = 2µ and h = 1. In both formulas ε1 should satisfy condition (42). For example, for µ = 1, the choice of ε1 = 10−10 guarantees accuracy ε ≤ 0.01 of the exact reach set approximation after 100 steps. Figure 2 (b) shows the regularized reach set at time step k = 10. Here we picked ε1 = 0.1, calculated the corresponding δ = 0.000025 and α = 0.45 and used the recurrence relation (43) to calculate the shape matrix of the external ellipsoid for different direction vectors l0 .
8
Internal ellipsoidal approximations
Consider the sum of k + 1 ellipsoids E(0, Hi ), 0 ≤ i ≤ k with positive semidefinite Hi . An internal ellipsoidal approximation of this sum is given by E(0, Q) ⊆
k
E(0, Hi ),
i=0
with Q = Q[k] =
k i=0
k 1/2 T
Si H i
i=0
1/2
Si H i
and any orthogonal matrices Si , 0 ≤ i ≤ k, SiT Si = I (see [1]). Q[k] satisfies the recurrence relation Q[k + 1] =
k+1
k+1 1/2 T 1/2 Si H i Si H i i=0 i=0
17
,
(49)
k
=
i=0
1/2
Si H i
1/2
+ Sk+1 Hk+1
= Q[k] + Hk+1 +
+
1/2
Sk+1 Hk+1
T
k
i=0 k i=0
k T i=0
1/2 T
Si H i
1/2
1/2
Si H i
1/2
+ Sk+1 Hk+1
1/2
Sk+1 Hk+1
Si H i
,
(50)
with Q[0] = H0 . The following tightness condition is true (see [1]). Theorem 8.1 Let H0 be positive definite. Then the inequality
l, Q[k]l =
l,
k
k 1/2 T 1/2 Si H i Si H i l i=0 i=0
≤
k
l, Hi l 1/2
2
,
(51)
i=0
which holds for any l ∈ Rn and any array of orthogonal matrices Si , 0 ≤ i ≤ k, turns into equality for a given l ∈ Rn iff for all Hi there exist numbers ηi such that 1/2
Si H i
1/2
l = ηi S0 H0 l,
0 ≤ i ≤ k.
(52)
Now let us return to system (1) with nonsingular matrices A[i], k0 ≤ i < k, and its reach set at time step k, X (k, k0 , E(x0 , X0 )) given by (6). As we can see, the reach set is a sum of k + 1 ellipsoids. If we define Q[k] as
1/2
Q[k] = Φ(k, k0 ) X0 S0T + ×
1/2
S0 X0
+
k−1
k−1 i=k0
(Φ(k0 , i + 1)R[i]ΦT (k0 , i + 1))1/2 SiT [k]
Si [k](Φ(k0 , i + 1)R[i]ΦT (k0 , i + 1))1/2 ΦT (k, k0 ),
(53)
i=k0
with Q[k0 ] = X0 , and in which S0 and Si [k], k0 ≤ i < k, are any orthogonal matrices, E(q[k], Q[k]) ⊆ E(Φ(k, k0 )x0 , Φ(k, k0 )X0 ΦT (k, k0 )) +E
k−1
Φ(k, i + 1)B[i]p[i],
i=k0
Φ(k, i + 1)R[i]ΦT (k, i + 1)
i=k0
= X (k, k0 , E(x0 , X0 )), with q[k] coming from (7). X (k, k0 , E(x0 , X0 )).
k−1
(54)
So, the ellipsoid E(q[k], Q[k]) lies inside the reach set
18
Definition 8.1 Ellipsoid E(q[k], Q[k]) ⊆ X (k, k0 , E(x0 , X0 )) is a tight internal approximation of the reach set X (k, k0 , E(x0 , X0 )) if there exists direction l such that ρ(±l | E(q[k], Q[k])) = ρ(±l | X (k, k0 , E(x0 , X0 ))). In order for E(q[k], Q[k]) to be a tight internal approximation of the reach set for a given direction l, we use Theorem 8.1 and come up with the following tightness condition that can be checked by direct calculation. Theorem 8.2 For a given time step k ≥ k0 and given direction l ρ(l | E(q[k], Q[k])) = ρ(l | X (k, k0 , E(x0 , X0 ))) iff S0 and Si [k], k0 ≤ i < k, satisfy the relation 1/2
Si [k](Φ(k0 , i + 1)R[i]ΦT (k0 , i + 1))1/2 ΦT (k, k0 )l = ηi [k]S0 X0 ΦT (k, k0 )l, with ηi [k] =
l, Φ(k, i + 1)R[i]ΦT (k, i + 1)l 1/2 . l, Φ(k, k0 )X0 ΦT (k, k0 )l 1/2
(55)
(56)
Let l[k] satisfy the good curve relation (17). In this case, expressions (55) and (56) simplify into 1/2 (57) Si (Φ(k0 , i + 1)R[i]ΦT (k0 , i + 1))1/2 l0 = ηi S0 X0 l0 , and ηi =
l0 , Φ(k0 , i + 1)R[i]ΦT (k0 , i + 1)l0 1/2 . l0 , X0 l0 1/2
(58)
Now ηi and Si , k0 ≤ i < k no longer depend on k. We are ready to write the recurrence relation for the shape matrix Q[k] of the internal ellipsoid E(q[k], Q[k] that touches the reach set X (k, k0 , E(x0 , X0 )) at points of the good trajectory specified by (17). Using definition (53) and relation (50) we get Q[k + 1] = A[k]Q[k]AT [k] + R[k]
+ Φ(k + 1, k0 ) Mk Sk (Φ(k0 , k + 1)R[k]ΦT (k0 , k + 1))1/2
+ (Φ(k0 , k + 1)R[k]ΦT (k0 , k + 1))1/2 SkT MkT ΦT (k + 1, k0 ),
(59)
with Q[k0 ] = X0 , and 1/2
Mk = X0 S0T +
k−1 i=k0
(Φ(k0 , i + 1)R[i]ΦT (k0 , i + 1))1/2 SiT .
(60)
The orthogonal matrices S0 , Sk , k ≥ k0 , are related by (57) and (58). The recurrence relation for the center of the internal approximating ellipsoid q[k] is the same as for external and is given by (22). 19
Theorem 8.3 Let l[k] satisfy the equation (17), then the ellipsoid E(q[k], Q[k]), where Q[k] and q[k] are the solutions of (59) and (22), is such that for all k ≥ k0 , the following is true: E(q[k], Q[k]) ⊆ X (k, k0 , E(x0 , X0 )) and ρ(l[k] | X (k, k0 , E(x0 , X0 ))) − ρ(l[k] | E(q[k], Q[k])) = 0. The control that drives the system along the good trajectory is indicated in (25) for those vectors l0 , for which l0 , Φ(k0 , i + 1)R[i]ΦT (k0 , i + 1)l0 = 0,
k0 ≤ i < k,
and the initial condition in (26). For every boundary point of the reach set X (k, k0 , E(x0 , X0 )) there exists direction l[k] and corresponding tight internal ellipsoid E(q[k], Q[k]). Thus, we arrive at the result similar to theorem 4.2. Theorem 8.4 The reach set X (k, k0 , E(x0 , X0 )) can be expressed as X (k, k0 , E(x0 , X0 )) =
E(q[k], Q[k]).
(61)
l[k]
Remark. The results of this section are also valid for singular matrices R[i], and no special regularization as in the case of external approximation is needed. Until now, in discussing internal ellipsoidal approximations, we assumed A[i] to be nonsingular. In case of singular A[i], the recurrence relations above lose their sense because matrix Φ(k0 , k) does not exist. It is still possible to build internal approximating ellipsoids using (49) even for the singular case, since the reach set of a linear discrete-time system for every time step k is nothing else but a finite sum of ellipsoids. However, the recurrence relations in the case of singular matrix A cease to exist, making the calculation of tight internal approximating ellipsoids a rather involved procedure because for every k the matrices Si need to be recomputed. Therefore, we suggest the following way out for singular case. Instead of the actual reach set X (k, k0 , E(x0 , X0 )) we consider the regularized reach set Xδ,α (k, k0 , E(x0 , X0 )) defined in (35) and constructed according to theorem 6.1. Matrices Aδ [i] are nonsingular, so are Rα [i], and the relations (57)-(60). Recall the example (48) from section 7. Figure 3 illustrates the internal ellipsoidal approximations √ ofthe regularized 1 3 1 2 √ 2 , and . reach set at k = 10 for three different initial directions, l0 = 3 0 − 12 2
20
1 0.8 0.6 0.4
x2
0.2 0 −0.2 −0.4 −0.6 −0.8 −1
−4
−3
−2
−1
0 x
1
2
3
4
1
Figure 3: Internal ellipsoidal approximations of the regularized reach set Xδ,α at k = 10.
9
Conclusion
This paper extends the results of [1] to linear discrete-time systems. We do not assume controllability of the system and also consider the situation when the state transition matrix Φ(k, k0 ) is singular. It was shown that in these situations it is possible to approximate the actual reach set with any given ε-accuracy by regularizing matrices B[i]P [i]B T [i] and A[i]. Theorem 6.1 provides the means of picking the appropriate regularization parameters. Internal ellipsoidal approximations work well even with singular matrices B[i]P [i]B T [i] and Φ(k, k0 ). However, in order to use the concept of ‘good curves’ and the recurrence relations (57)-(60), we still need to perform the regularization.
References [1] A. Kurzhanski and P. Varaiya, “On ellipsoidal techniques for reachability analysis,” Optimization Methods and Software, vol. 17, pp. 177–237, 2000. [2] N. Giorgetti and G. Pappas, “Bounded model checking for hybrid dynamical systems,” in IEEE Conference on Decision and Control, Seville, Spain, 2005, submitted. [3] I. Mitchell and C. Tomlin, “Level set methods for computation in hybrid systems,” in Hybrid Systems: Computation and Control, ser. Lecture Notes in Computer Science, N. Lynch and B. Krogh, Eds. Springer, 2000, vol. 1790, pp. 21–31. [4] “Level set toolbox,” www.cs.ubc.ca/˜mitchell/ToolboxLS.
21
[5] M. Kvasnica, P. Grieder, M. Baoti´c, and M. Morari, “Multi parametric toolbox (mpt),” in Hybrid Systems: Computation and Control, R. Alur and G. Pappas, Eds. Springer, 2004, vol. 2993, pp. 448–462. [6] “Mpt homepage,” control.ee.ethz.ch/˜mpt. [7] T. Motzkin, H. Raiffa, G. Thompson, and R. Thrall, “The double description method,” in Conttributions to Theory of Games, H.W.Kuhn and A.W.Tucker, Eds. Princeton University Press, 1953, vol. 2. [8] “Cdd/cdd+ homepage,” www.cs.mcgill.ca/˜fukuda/soft/cdd home/cdd.html. [9] “d/dt homepage,” www-verimag.imag.fr/˜tdang/ddt.html. [10] E. Asarin, O. Bournez, T. Dang, and O. Maler, “Approximate reachability analysis of piecewise linear dynamical systems,” in Hybrid Systems: Computation and Control, ser. Lecture Notes in Computer Science, N.Lynch and B.H.Krogh, Eds. Springer, 2000, vol. 1790, pp. 21–31. [11] D. Avis, D. Bremner, and R. Seidel, “How good are convex hull algorithms?” Computational Geometry: Theory and Applications, vol. 7, pp. 265–301, 1997. [12] A. Girard, “Reachability of uncertain linear systems using zonotopes,” in Hybrid Systems: Computation and Control, ser. Lecture Notes in Computer Science. Springer, 2005, to appear. [13] “Zonotope methods on wolfgang k¨ uhn homepage,” www.decatur.de. [14] “CheckMate homepage,” www.ece.cmu.edu/˜webk/checkmate. [15] O. Stursberg and B. Krogh, “Efficient representation and computation of reachable sets for hybrid systems,” in Hybrid Systems: Computation and Control, ser. Lecture Notes in Computer Science, O.Maler and A.Pnueli, Eds. Springer, 2003, vol. 2623, pp. 482–497. [16] “Requiem homepage,” www.seas.upenn.edu/˜hybrid/requiem/requiem.html. [17] G. Lafferriere, G. Pappas, and S. Yovine, “Symbolic reachability computation for families of linear vector fields,” Journal of Symbolic Computation, vol. 32, pp. 231–253, 2001. [18] E. Kostousova, “Control synthesis via parallelotopes: optimization and parallel computations,” Optimization Methods and Software, vol. 14, no. 4, pp. 267–310, 2001. [19] P. Varaiya, “Reach set computation using optimal control,” Proc. of KIT Workshop on Verification on Hybrid Systems. Verimag, Grenoble., 1998. [20] A. Kurzhanski and P. Varaiya, “Dynamic optimization for reachability related problems,” Journal of Optimization, Theory and Applications, vol. 108, no. 2, pp. 227–251, 2001.
22
