The paper suggests the modification of the generalized recursion algorithm of the exact reachable sets computation for the linear discrete dynamic 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.
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)
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: 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:
x(0) = x0 ∈ X(0) ⊂ Rn, u(t) ∈ P(t) ⊂ Rp ∀t ∈ 0, T − 1, 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
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).
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 (t, G+(t)) ∈ 0, T − 1 × 2Rn , 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
Gb x(t + 1) = {xˆ(t + 1) ∈ Rn | xˆ(t + 1) = A(t)x(t), x(t) ∈ Γn (G+(t))} ,
Gb u(t + 1) = {yˆ(t + 1) ∈ Rn | yˆ(t + 1) = B(t)u(t), u(t) ∈ Γp (P(t))} ,
Gb +(t + 1) = nvˆ(t + 1) ∈ Rn | vˆ(t + 1) = xˆ(t + 1) + yˆ(t + 1), xˆ(t + 1) ∈ Gb x(t + 1), yˆ(t + 1) ∈ Gb u(t + 1)o , where Gb +(t + 1) is the Minkowski sum of the sets Gb x(t + 1) and Gb 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 Gb +(t + 1) = vˆ(i)(t + 1) i∈1,m ⊂ G (t, G+(t); t + 1) we define the set of all its vertices Γn
Gb +(t + 1) = vˆ(i)(t + 1) i∈1,k , k ≤ m.
Taking into account that, according to [Shorikov97, Tyulyukin93], the following statement holds true: conv Gb +(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. (2) (3) (4) In the works [Shorikov97, Tyulyukin93] it was shown that the set Gb +(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: f =
X λj → min,
j s.t.
X λj vˆ(j)(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 there exists feasible basic solution (FBS), then the point corresponding the vector vˆ(i)(t + 1) is not the vertex of the polytope conv Gb +(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
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 = Pkn=1a1na1k1 + 1 ··· ··· ··· Pkn=1b1nbk + 1 ··· ··· ··· Pkn=1ana1mkm + 1, bk ≥ 0, k ∈ 1, n.
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 could be presented in the form: b0 1 . . . b0 n b0n+1 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
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 Ox1x2 (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:
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.
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
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). 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).
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.
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 algorithm [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
The research was supported by Russian Foundation for Basic Research (Project 15-01-02368).