=Paper=
{{Paper
|id=Vol-1825/p11
|storemode=property
|title= A modification of the generalized recursion method of the linear control systems reachable sets computation
|pdfUrl=https://ceur-ws.org/Vol-1825/p11.pdf
|volume=Vol-1825
|authors=Andrew F. Shorikov,Vladimir V. Bulaev
}}
== A modification of the generalized recursion method of the linear control systems reachable sets computation ==
https://ceur-ws.org/Vol-1825/p11.pdf
A modification of the generalized recursion method of
the linear control systems reachable sets computation
Andrei F. Shorikov1 Vladimir V. Bulaev2
afshorikov@mail.ru bulaev1991@mail.ru
1–Ural Federal University (Yekaterinburg, Russia)
2–JSC ”Scientific and Production Association of Automatics” (Yekaterinburg, Russia)
Abstract
The paper suggests the modification of the generalized recursion algo-
rithm of the exact reachable sets computation for the linear discrete dy-
namic systems. The examples of the finite convex hull search problem
and the reachable sets computation problem demonstrate comparative
analysis of the origin and modified methods.
Keywords: reachable set, discrete-time dynamic system, simplex
method, convex hull.
1 Introduction
In the modern control theory there are lots of problems which require computation and analysis of the control
system reachable sets [Chernousko94, Krasovskii68, Krasovskii88]. In principle, the methods of the reachable
sets computation could be separated into two approaches. The first approach involves computation of the exact
reachable sets, comprising only those dynamic systems states in which it could be transformed for the finite period
of time. These methods were developed in the research works [Krasovskii88, Lasserre91, Lotov72, Shorikov97,
Tyulyukin93]. The second approach includes some approximation of the reachable set, for example, with the
use of ellipsoids [Chernousko94, Kurzhanski02] or parallelotopes [Kostousova01]. Approximation methods up to
date are paid a lot of attention [Asarin00, Girard05, Kurzhanskiy11], since they allow decreasing computation
costs. However approximation methods give only unfaithful representation and in some cases rough idea about
reachable sets.
This paper suggests the modification of the generalized recursion algorithm of the exact reachable sets com-
putation for the linear discrete-time dynamic systems developed in the works [Shorikov97, Tyulyukin93].
2 Dynamic system description
Consider on the integer-value period of time t ∈ 0, T = {0, 1, . . . , T }, T ∈ N (N is the set of all natural numbers)
linear control system with the dynamics described by the discrete vector-matrix recursion equation of the form:
x(t + 1) = A(t)x(t) + B(t)u(t), t ∈ 0, T − 1, (1)
where x(t) is the state vector (phase vector), x(t) ∈ Rn (from now on Rn is the n-dimensional Euclidean space of
the column-vectors); u(t) is the control vector, u(t) ∈ Rp ; A(t) is the system state matrix, A(t) ∈ Rn×n ; B(t) is
the control matrix, B(t) ∈ Rn×p .
Copyright c by the paper’s authors. Copying permitted for private and academic purposes.
In: G.A. Timofeeva, A.V. Martynenko (eds.): Proceedings of 3rd Russian Conference ”Mathematical Modeling and Information
Technologies” (MMIT 2016), Yekaterinburg, Russia, 16-Nov-2016, published at http://ceur-ws.org
88
For the description of the considered dynamic systems class we introduce the assumptions about the initial
state x(0) of the system (1) and control vector u(t).
Assumption 2.1 Phase vector initial states of the system (1) satisfy the defined geometric constraint:
x(0) = x0 ∈ X(0) ⊂ Rn , (2)
where X(0) is the convex, closed and limited polytope with the finite number of vertices.
Assumption 2.2 Control vector values satisfy the defined geometric constraint:
u(t) ∈ P(t) ⊂ Rp ∀t ∈ 0, T − 1, (3)
where P(t) is the convex, closed and limited polytope with the finite number of vertices.
Under the conditions of the faithful Assumptions 2.1 and 2.2 we introduce the reachable set definition for the
discrete-time control system (1) – (3) [Krasovskii68, Shorikov97].
Defenition 1 The reachable set of the control system phase states (1) – (3) on the final moment of time T ,
which is correspondent to the pair (0, X(0)), is the set defined as following:
G (0, X(0); T ) = x(T ) | x(T ) ∈ Rn , x(t + 1) = A(t)x(t) + B(t)u(t), x(0) ∈ X(0), u(t) ∈ P(t), t ∈ 0, T − 1 .
3 Generalized recursion method of the reachable sets computation
In the works [Shorikov97, Tyulyukin93] it was shown that the reachable set corresponding to the Definition 1
represents the convex, closed, limited polytope with the finite number of vertices in the space Rn . Besides, for
such reachable set the following recursion (semigroup) property holds true:
G (0, X(0); T ) = G (t, G+ (t); T ) ∀t ∈ 1, T − 1,
where G+ (t) = G (0, X(0); t) is the reachable set for the moment of time t, corresponding to the pair (0, X(0))
which is the convex, closed, limited polytope with the finite number of vertices in Rn .
In this case, the set G (0, X(0); T ) search problem reduces to the recursion computation sequence of the
one-step reachable sets:
G (t, G+ (t); t + 1) , t ∈ 1, T − 1, G+ (0) = X(0). (4)
Consequently, the basic auxiliary problem in defining the set G (0, X(0); T ) is the problem of the reachable sets
computation (4), for instance, with the use of its vertices sets determination.
Hence we define the approach towards the defined basic auxiliary problem solution, relying on the generalized
recursion algorithm [Shorikov97].
The set of the feasible phase states x(t + 1) of the considered control system (1) – (3), comprising all the
reachable set G (t, G+ (t); t + 1) vertices on the moment of time (t + 1), corresponding to the pair
n
(t, G+ (t)) ∈ 0, T − 1 × 2R , is defined according to the following algorithm.
Step 1. Forming of the set Γn (G+ (t)) of all polytope G+ (t) vertices.
Step 2. Forming of the set Γp (P(t)) of all polytope P(t) vertices.
Step 3. Defining the following sets, characterizing free and disturbed motion
b x (t + 1) = {x̂(t + 1) ∈ Rn | x̂(t + 1) = A(t)x(t), x(t) ∈ Γn (G+ (t))} ,
G
Gb u (t + 1) = {ŷ(t + 1) ∈ Rn | ŷ(t + 1) = B(t)u(t), u(t) ∈ Γp (P(t))} ,
n o
b + (t + 1) = v̂(t + 1) ∈ Rn | v̂(t + 1) = x̂(t + 1) + ŷ(t + 1), x̂(t + 1) ∈ G
G b x (t + 1), ŷ(t + 1) ∈ G
b u (t + 1) ,
where G b + (t + 1) is the Minkowski sum of the sets G b x (t + 1) and G
b u (t + 1), being the set of the reachable set
points, among which are both internal and frontier points.
Step 4. For the defined set G b + (t + 1) = v̂ (i) (t + 1) ⊂ G (t, G+ (t); t + 1) we define the set of all its
i∈1,m
vertices Γn G b + (t + 1) = v̂ (i) (t + 1) , k ≤ m.
i∈1,k
Taking into account that, according to [Shorikov97, Tyulyukin93], the following statement holds true:
conv G b + (t + 1) = G (t, G+ (t); t + 1) = G+ (t + 1),
then the vertices of the set Gb + (t + 1) are defining the control system (1) – (3) reachable set for the moment
(t + 1), which is the convex, closed, limited polytope in Rn .
89
3.1 Vertices search problem
In the works [Shorikov97, Tyulyukin93] it was shown that the set G
b + (t + 1) vertices search problem is reduced
to the solution of m linear mathematical programming (LP) problems. In accordance with this fact, we are
formulating the following optimization problem.
Problem 1 (LP1) For the fixed i ∈ 1, m and the set of variables λj , j ∈ 1, m, j 6= i, λi ∈ R1 , it is required to
solve the following LP problem:
X
f= λj → min,
j
X
(j)
s.t. λj v̂ (t + 1) = v̂ (i) (t + 1),
j
X
λj = 1, λj ≥ 0.
j
It should be noted that the cost function f in the LP1, in general terms, could have any form, since in the
process of this problem solving we are interested only in one question, if the considered problem has the basis
feasible solution or not. The constraints system properties of the LP1 follow that if it is consistent, namely
(i)
there exists feasible basic
solution
(FBS), then the point corresponding the vector v̂ (t + 1) is not the vertex
of the polytope conv G b + (t + 1) , since it could be presented as the convex combination of the other points.
Otherwise, the considered point is the reachable set G+ (t + 1) vertex.
As a result, solving m LP1 problems, we could find the whole set of vertices Γn (G+ (t + 1)) which does present
the set of all the reachable set G (t, G+ (t); t + 1) vertices for the dynamic system (1) – (3) for the moment of
time (t + 1).
In the original version of algorithm [Shorikov97], for the solution of LP1 it is suggested to use modified simplex
method (MSM) [Papadimitriou82, Yudin69].
4 Modified method of the reachable sets computation
In terms of computational complexity, in the considered algorithm of reachable set computation significant role is
played by the MSM, since it has exponential complexity in relation to the dimension of constraints set – (n + 1)
[Papadimitriou82]. Consequently, the efficiency of the whole algorithm, in general, is directly defined by the
efficiency of the LP1 solution algorithm.
The modification idea of original reachable set computation method [Shorikov97, Tyulyukin93] is in reducing
the number of solved LP1 problems due to the more complete use of the geometric properties of simplex method.
Consider the process of the LP1 problem solution. In order to find the FBS of the LP constraints system we
use the artificial basis technique [Yudin69]. The initial simplex table for some verifiable point v̂ (i) (t + 1), i ∈ 1, m
looks as following:
a11 ··· b1 ··· a1m
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
M =
an1 ··· bn ··· anm , bk ≥ 0, k ∈ 1, n.
Pn 1 ··· P 1 ··· P 1
n n
k=1 ak1 + 1 · · · k=1 bk + 1 · · · k=1 akm + 1
Where vector b = {bk }k∈1,n corresponds to the coordinates of the verified point v̂ (i) (t + 1) and vectors
a(j) = {akj }k∈1,n , j ∈ 1, m, j 6= i, correspond to the other points. The next to last row of the initial matrix
M represents the constraint, related to the definition of the convex combination. The last row defines the
substitution cost of the corresponding simplex table columns.
Assume that the considered LP problem has the FBS, namely the verified point v̂ (i) (t + 1) is internal. As
the choice criteria for the marker column in the simplex table M we use the method of non-basis gradient
[Papadimitriou82].
As a result, at the final stage of the simplex method first phase, without loss of generality, the matrix M
90
could be presented in the form:
0
b01 a01l
a11 ··· ··· 1 ··· 0 0
.. .. .. .. . . .. ..
. ··· . ··· . . . . . 0
M = a , b ≥ 0, k ∈ 1, n + 1, m − (n + 1) ≤ l < m.
0
0 n1 ··· b0n ··· a0nl 0 ··· 1 0 k
an+1,1 · · · b0n+1 ··· 0
an+1,l 0 ··· 0 1
0 ··· 0 ··· 0 0 ··· 0 0
Define B = {a(j) }j∈l+1,m as the set of basis vectors. In this case, vector b could be presented as the convex
combination of the basis vectors a(i) ∈ B, i ∈ 1, m − l. Besides, among other points, not corresponding to
the basis vectors, could be such points which could be expressed as the convex combination of basis vectors
a(i) . Namely, according to [Papadimitriou82], we should find the columns of the matrix M with the following
properties:
n+1
X
a0kj = 1, a0kj ≥ 0 ∀k ∈ 1, n + 1, a(j) ∈
/ B, j ∈ 1, l.
k=1
Geometrically, this means that these points could be comprised by simplex with the dimension not more
than (n + 1), with vertices relying on the basis vectors a(i) ∈ B. Consequently, they could be escaped from the
following consideration, reducing the number of the solved LP1 problems. This is the first part of modification
of the original method [Shorikov97].
The second part of modification claims that if the verified point v̂ (i) (t + 1) is internal, then the points,
corresponding to the set of basis vectors a(i) ∈ B, represent the subset of the all reachable set vertices. However
in this part of modification there is some deviation from the classic method of non-basis gradient for the choice
of marker columns direction in the simplex table M .
The following is the example for the two-dimensional subspace case.
4.1 Example
Consider the set of points A(1, 1), B(2, 2), C(0, 3/2), D(1, 1/2), E(3, 1), F (1, 3), G(−1, 0) on the plane Ox1 x2
(Figure 1). It is required to find a convex hull for the defined set of points. Let A(1, 1) be the verified point.
Then the initial simplex table has the form:
1 2 0 1 3 1 −1
1 2 3/2 1/2 1 3 0
M = 1 1
.
1 1 1 1 1
3 5 5/2 5/2 5 5 0
Figure 1: Illustration for the Example
If one implements the search of the FBS for such problem with the use of the classical method of non-basis
gradient, namely if at each iteration from the left to the right one chooses columns with the maximal substitution
cost (values in the last row), then as a result the basis vectors occur to be corresponding to the points B(2, 2),
C(0, 3/2), D(1, 1/2). Namely point A(1, 1) will be comprised by the triangle BCD, vertices of which are not the
vertices of analyzed set of points.
91
Note that the equality of the substitution costs of the second, fifth and sixth columns in initial simplex table
M is explained by the fact that corresponding points are situated on the same line x1 + x2 = 4, characterizing
the level of the maximal substitution costs. After the choice of the first basis vector, corresponding to the point
B(2, 2), and implementation of the first iteration, the equality of the substitution costs holds for the points
C(0, 3/2), F (1, 3), G(−1, 0) which are situated on the line −3x1 + 2x2 = 3. These lines we will called specific
lines. The equation of the second specific line could be derived from the equality condition for the substitution
costs after the first iteration and the condition x1 + x2 ≤ 4. This follows that equation of specific line (except for
the first one) is defined by the basis vector choice. After the choice of the second basis vector, corresponding to
the point C(0, 3/2), the equality of the substitution costs holds for the points D(1, 1/2), E(3, 1), G(−1, 0) which
are situated on the line −x1 + 4x2 = 1. The equation of this specific line is derived from the equality condition
for the substitution costs after the second iteration and the conditions x1 + x2 ≤ 4 and −3x1 + 2x2 ≤ 3.
In order to avoid such situations, when the verified point could be comprised by the simplex, relying on the
points which are not vertices, it is enough to choose from the columns with the equal substitution costs the
marker column with the maximal norm function. If this property is taken into account, then the basis vectors
occur to be the vectors, corresponding to the vertices E(3, 1), F (1, 3), G(−1, 0). Namely, the verified point
A(1, 1) will be comprised by the triangle EF G. Besides, from the consideration the points B(2, 2), C(0, 3/2),
D(1, 1/2) could be escaped, since, as the verified point, they will be the convex combinations of these vertices.
As a result, the convex hull search problem is solved for the one MSM procedure, including only three iterations.
In this way, in the example we have clarified the sense of the second part of the original modification method
[Shorikov97]. Analogously to the given example, for the higher dimensions spaces the concept of ”the specific
line” is expanded to the concepts of ”the specific plane” or ”the specific hyperplane”. However, it should be
noted that such specific cases are quite rare at the practice, since they require specific conditions.
It should be mentioned, that the described modifications are applied only in the cases when the verified point
is internal. In the case when it is vertex, it is necessary to include the verified point into the list of vertices and
turn to the verification of the next point.
5 Comparative analysis of the original method and its modification
In order to analyze the improvement of original method efficiency [Shorikov97] due to the implementation of the
described modifications, there is accomplished the comparative analysis, consisting of the two stages. At the first
stage the evaluation was performed at the example of finite random convex hull search problem. At the second
stage the speed of work was estimated at the example of the reachable sets (4) computation for the particular
models of the discrete-time control systems (1) – (3). These are the indicators of the speed estimates: computation
time and the total number of the MSM iterations.
Turn to the first stage of the comparative analysis. For the productivity estimation the problems of convex
hull computation were solved for three typical sets of random points (Figure 2).
Figure 2: Convex hulls of the typical sets in R2
1. Type Normal is the random set of points with the normal law of distribution (expected mean M = 0,
standard deviation σ = 0.5).
2. Type Uniform is the random set of points with the uniform law of distribution (M = 0, σ = 1).
3. Type Sphere is the random set of points, normally distributed relative to the sphere surface with the radius
of R = 5 (M = 4.75, σ = 0.2).
92
Figure 3: Results of the numerical experiments
On the Figure 3 dotted graph means results obtained with the use of the original algorithm, solid graph
denotes results obtained with the use of the modified algorithm. The first graph reflects the computing time of
algorithms depending on the Euclidean space dimension n. The second graph shows the total MSM iterations
number growth in the logarithmic scale.
Table 1: The results of the reachable sets calculation
Calculation time, sec Total number of MSM iteration
Time moment Original Modified Original Modified
1 0,002865 0,003340 46 46
2 0,008091 0,007398 211 192
3 0,019112 0,015168 609 500
4 0,041457 0,035045 1 468 1 198
5 0,093163 0,076567 3 036 2 454
6 0,171463 0,119941 5 631 4 399
7 0,267551 0,227781 9 675 7 487
8 0,450815 0,267551 15 300 11 740
9 0,646412 0,448515 23 136 17 543
10 0,959281 0,657190 33 990 25 600
11 1,424411 0,958765 48 266 36 177
12 2,104862 1,433354 67 092 50 086
13 2,858660 1,930592 91 507 68 005
14 3,804880 2,573760 121 407 89 804
15 5,466993 3,583817 158 399 116 497
In order to demonstrate the second stage of the comparative analysis we show the results of the reachable sets
computation for the following model of the discrete dynamic fourth order model:
1 0 0.03 0 −0.0024
0.014 1 0.0001 0.03 x(t) + 0.0051 u(t), x(t) ∈ R4 , u(t) ∈ R1 , t ∈ 0, 14
x(t + 1) =
0.114 0 1 −0.0002 −0.157
0.915 0 0.014 1 0.339
|x1 (0)| ≤ 0.052, |x2 (0)| ≤ 5, |x3 (0)| ≤ 0.018, |x4 (0)| ≤ 1, |u(t)| ≤ 0.14.
The results of the comparative analysis second stage for this discrete-time dynamic system are presented in
93
the Table 1.
6 Conclusion
Under the results of the numerical experiments, presented at the Figure 3, it could be concluded that modified
method in the frames of the convex hull computation problem works significantly faster than the original algo-
rithm [Shorikov97]. It is especially seen for the low dimensional spaces (n = 2, . . . , 5). The convergence of the
speed for two considered algorithms with higher n is explained by the fact that they involve the MSM, which is
known [Papadimitriou82] to have exponential complexity relative to n.
The results of the reachable sets computation problem (presented in Table 1) show that the modification of
the generalized recursion method allow significantly reducing the number of MSM iterations and, consequently,
significantly reducing the time of computations.
6.0.1 Acknowledgements
The research was supported by Russian Foundation for Basic Research (Project 15-01-02368).
References
[Asarin00] E. Asarin, O. Bournez. Approximate reachability analysis of piecewise linear dynamical systems.
Hybrid Systems: Computation and Control, 1790:21–31, 2000.
[Chernousko94] F. L. Chernousko. State estimation for dynamic systems, CRC Press, Bocaraton, Florida, 1994.
[Girard05] A. Girard. Reachability of uncertain linear systems using zonotopes. Hybrid Systems: Computation
and Control, 3414:291–305, 2005.
[Kostousova01] E. K. Kostousova Control synthesis via parallelotopes: optimization and parallel computations.
Optimization Methods and Software, 14:267–310, 2001.
[Krasovskii68] N. N. Krasovskii. Theory of control of motion: Linear systems (Russian). Moscow, Nauka, 1968.
[Krasovskii88] N. N. Krasovskii, A. I. Subbotin. Game-Theoretical Control Problems. New York, Springer-Verlag,
1988.
[Kurzhanski02] A. B. Kurzhanski, P. Varaiya. On ellipsoidal techniques for reachability analysis. Optimization
Methods and Software, 17:177–207, 2002.
[Kurzhanskiy11] A. A. Kurzhanskiy, P. Varaiya. Reach set computation and control synthesis for discrete-time
dynamical systems with disturbances. Automatica, 47:1414–1426, 2011.
[Lasserre91] J. B. Lasserre. Reachable and controllable sets for two-dimensional, linear, discrete-time systems.
Journal of Optimization Theory and Applications, 70:583–595, 1991.
[Lotov72] A. V. Lotov. Numerical method of constructing attainability sets for a linear control system (Russian).
Computational Mathematics and Mathematical Physics, 12(3):785–788, 1972.
[Papadimitriou82] C. H. Papadimitriou, K. Steiglitz Combinatorial optimization: algorithms and complexity.
Prentice-Hall, Englewood Cliffs, New Jersey, 1982.
[Shorikov97] A. F. Shorikov. Minimax estimation and control in discrete-time dynamic systems (Russian).
Yekaterinburg, Ural State University, 1997.
[Tyulyukin93] V. A. Tyulyukin, A. F. Shorikov. The solution algorithm of terminal control problem for linear
discrete-time system (Russian). Automatics and Telemechanics, 4:115–127, 1993.
[Yudin69] D. B. Yudin, E. G. Golshtain. Linear Programming: theory, methods and approaches (Russian).
Moskow, Nauka, 1969.
94