Reachable set computation for solving fuel consumption terminal control problem using zonotopes Andrey F. Shorikov1 Vitaly I. Kalev2 afshorikov@mail.ru butahlecoq@gmail.com 1–Ural Federal University (Yekaterinburg, Russia) 2–JSC ”Scientific and Production Association of Automatics” (Yekaterinburg, Russia) Abstract In this paper, fuel consumption terminal control problem of liquid- propellant rocket is considered as an optimal terminal control problem of linear discrete-time dynamical system. Control of the system is sup- posed bounded. Also, we assume there are no external disturbances affect on the system. For solving optimal terminal control problem we introduce the problem of exact reachable set computation. Usu- ally for discrete-time dynamical systems the exact reachable sets have polytope representation. This paper provides an algorithm of exact reachable set computation using zonotopes to represent the exact reach- able sets of discrete-time dynamical system. Zonotopes are the special class of polytopes, which are usually used for approximations of reach- able sets. Zonotopes have an efficient property: they are closed under linear transformations and Minkowski sum. Performance of proposed algorithm is shown in comparison to exact reachable set computation algorithm using polytopes at the modeling example of launch vehicle two-tank sixteen-step discrete-time dynamical system in R4 . Keywords: Reachable Set, Polytopes, Zonotopes, Terminal Control, Optimal Control, Launch Vehicle, Fuel Consumption Control. 1 Introduction Most of the modern space missilery plants which have liquid-propellant engine (e.g. launch vehicles, long-range missiles, spacecrafts) involve four basic control problems: guidance, navigation, control (stabilization) and fuel consumption control. The fuel consumption control problem of liquid-propellant vehicle is to synchronously consume fuel components in two tanks: oxidizer tank and fuel tank. Furthermore, it is needed to completely consume fuel supply at a given time. Thus, the fuel consumption control problem could be considered as an optimal terminal control problem, where the quality criterion is to minimize deviation of fuel rate and amount of fuel from their desired values at final time instant. Over the past 40 years a large number of papers have been devoted to solving fuel consumption terminal control problem of liquid-propellant vehicle. The results of the long investigation of this problem are summarized in [AI13]. Copyright © 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 103 At the present time the optimization problems of dynamical system terminal control with known probabilistic characteristics of a-priori uncertain parameters are investigated well enough. However, in real life, there is an absence of such probabilistic characteristics of a-priori unknown initial state vector [Sh93]. To handle the optimal control of complex plants with terminal cost function one needs to take it into account. In order to address this problem we need to study system dynamics in the time domain. The central concept arising in such studies is the reachable set – such set of states to which the dynamical system can be steered using feasible controls. This paper is a continuation of our research (see [SK16-1], [SK16-2]) that aimed to solve the terminal control problem of fuel consumption of the liquid-propellant launch vehicle. The previous papers [SK16-1], [SK16-2] are devoted to forming of linear discrete-time model of launch vehicle dynamics in terms of propulsion system dynamics and to modeling of obtained linear discrete-time system on PC. In this paper, we present algorithm for exact reachable set computation that allow us to solve above optimal terminal control problem later. In order to describe exact reachable set we use two representations: V-polytope (vertex representation of polytope) (e.g. see [Zi98], [Sh93], [Pa95]) and zonotope (see [Fu04], [Co03]). The proposed reachable set computation algorithm based on ideas stated in [Sh93] and has one difference: to represent exact reachable sets of linear discrete-time dynamical system new algorithm uses zonotopes instead of V-polytopes. The paper is organized as follows. First, we briefly describe the fuel consumption control problem. Then, we introduce the mathematical notion of zonotope. Afterwards, we provide algorithm for exact reachable set computation using zonotopes. Finally, performance of proposed algorithm is shown in comparison to exact reachable set computation algorithm [Sh93] using polytopes at the modeling example of launch vehicle two-tank sixteen-step discrete-time dynamical system in R4 . 2 Fuel consumption terminal control problem Above we said that on-board fuel consumption control system, implemented on some liquid-propellant launch vehicles, is required for increasing its energy potential. To achieve this purpose it is necessary to address the main problem that is completely and synchronously consume fuel components at a given time. To solve above the main problem we present the model of linear discrete-time dynamical system, constraints of state and control vector, and quality criterion. 2.1 Mathematical model At the given time interval 0, T = {0, 1, . . . , T } we consider linear discrete-time system that represent the dynamics of liquid-propellant launch vehicle propulsion system with two fuel tanks shown on figure 2.1. The system has the following state-space representation: x(t + 1) = A(t)x(t) + B(t)u(t), (1) where x(·) ∈ Rn is the state vector; u(·) ∈ Rr is the control (input) vector; A(·) ∈ Rn×n is the state matrix; B(t) ∈ Rn×r is the input matrix; t is the discrete time, t ∈ 0, T − 1. The state vector of this dynamical system includes masses and rates of fuel as the state coordinates. Detailed procedure of obtaining such linear discrete-time model of propulsion system and fuel tanks of launch vehicle described in [SK16-1]. 2.2 System constraints The initial state vector of the plant is determined by variety of propulsion system initial setting-up errors (e.g. filling error, throttle error, combustion chamber pressure error, etc). Since probabilistic characteristics of propulsion system initial setting-up errors are unknown, but only their upper bounds, one can consider the following initial set of states: x(0) ∈ X0 ⊂ Rn . (2) The requirement of propulsion system safe work (explosion prevention) constrains the control vector by the set: u(t) ∈ P(t) ⊂ Rp . (3) 104 Figure 1: Fuel consumption control system of launch vehicle‘s first stage 2.3 Terminal quality criterion The main problem is to completely and synchronously consume fuel components at a given time T . Moreover, it should be taken into account that the better values of fuel rates in propulsion system at final time T is to be equal to nominal (effective) values. Therefore, it is necessary to minimize deviation of state vector x(T ) at time T (fuel rate and amount of fuel) from desired values xd , i.e. to minimize the following Euclidean norm Φ : Rn × Rr×(T −1) → R1 : Φ = σ x(T ) − xd  , (4) n   where x(T ) = x x(0), {u(t)}t∈1,T −1 ∈ Rn – final state vector of plant (1); xd – desired final state vector of plant (1); σ – scale coefficient vector; || · || – Euclidean norm in Rn . Accordingly, the optimal terminal control problem of launch vehicle fuel consumption leads to finding such control sequence u(e) (·) = u(e) (t) t∈1,T −1 ∈ Rr×T −1 , ∀t ∈ 1, T − 1 : u(e) (t) ∈ P(t),  that steers the system (1) at time T from initial set (2) to such final set, in which functional (4) possesses the minimum value. 3 Zonotopes: definition and properties Zonotopes are a special class of convex polytopes. Zonotope is usually defined as the linear image of a p- dimensional hypercube in Rn . Equivalently, a zonotope is a Minkowski sum of a finite set of line segments. In this paper, we keep the following definition: Definition 1 (Zonotope). A zonotope Z is a set such that: i=p ( ) X n Z= x∈R :x=c+ xi ri , −1 6 xi 6 1 , i=1  where c, r1 , . . . , rp are vectors in Rn . We note Z = c, hr1 , . . . , rp i .  Zonotope Z = c, hr1 , . . . , rp i is always centrally symmetric and the point c ∈ Rn is the center of Z. The set of vectors r1 , . . . , rp is called the set of generators of zonotope Z. On figure 2 we show a zonotope construction with three generators in R2 . 105 Figure 2: Construction of zonotope with three generators in R2 Zonotopes have long been studied in combinatorial geometry [Zi98]. At the present time zonotopes are usually applied to approximate the reachable sets of uncertain dynamical systems (e.g. [Co03], [Gi05], [ASB10]). In this paper, zonotopes are used to compute exact reachable set of linear discrete-time dynamical system (1). Finally, let us show two main properties of zonotopes:  1. Zonotopes are closed under linear transformation. Let L be a linear map and Z = c, hr1 , . . . , rp i a zonotope, i=p ( ) X  LZ = Lx : x = c + xi ri , −1 6 xi 6 1 = Lc, hLr1 , . . . , Lrp i , p ∈ N. i=1 One can compute the image of a zonotope by linear map in linear time with regard to the order of the zonotope. Order of the zonotope is p/n, where p – number of generators and n – space dimension.   2. Zonotopes are closed under Minkowski sum. Let Z1 = c1 , hr1 , . . . , rp i and Z2 = c2 , hl1 , . . . , lq i be two zonotopes,  Z1 + Z2 = c1 + c2 , hr1 , . . . , rp , l1 , . . . , lq i , p ∈ N, q ∈ N. Hence, one can compute the Minkowski sum of two zonotopes by simple concatenation of two lists. 4 Reachable sets: definition and properties Before we present an algorithm of a reachable set computation, let us introduce the definition and some properties of a reachable set of system (1). Following definition of a reachable set is based on [Sh93].  Definition 2 (Reachable set). Reachable set of system (1)–(3), corresponding to pair i, X(i) ∈ 0, T − 1 × n ×2R , at time j (j ∈ i + 1, T ) is a set G i, X(i), j = x(j) : x(j) ∈ Rn ∀t ∈ i, j − 1,   x(t + 1) = Ax(t) + Bu(t), x(i) ∈ X(i), u(t) ∈ P(t) . In [Sh93] it is assumed that sets X0 and P(t) are the convex polytopes. Taking into account this assumption and considering the linearity of system (1), we can present the following properties of reachable set.  Property 1 (Compactness and Convexness). A set G 0, X0 , t for all t ∈ 1, T is a convex polytope in Rn . Property 2 (Evolutionary property).   G 0, X0 , t + 1 = G t, X(t), t + 1 , where X(t) = G 0, X0 , t), t ∈ 1, T − 1.  It is important to note that every reachable set G 0, X0 , t , ∀t ∈ 1, T is exact reachable set, since it includes only such states of system in which system can be steered by feasible control u(t) ∈ P(t). From the evolutionary property we have the following.  Corollary 1. Multi-step problem of reachable set G 0, X0 , T computation leads to solving of recurrent sequence  of one-step problems of reachable set X(t + 1) = G t, X(t), t + 1 , t ∈ 0, T − 1 computation. 106  Also, in [Sh93] it is shown that reachable sets G 0, X0 , t , ∀t ∈ 1, T can be represented as the set of all vertices of the polytope, that is as so-called V-polytope (see also [Zi98], [Fu04]). In this paper, to represent reachable set of system (1)–(3) we use zonotopes. Since the problem is to compute exact reachable set of system (1), it is necessary to do the following assumptions. Assumption 1. The initial set of states X0 is zonotope, i.e. it can be represented as  X0 = c1 , hr1 , . . . , rp i , p ∈ N. Assumption 2. The constraint of control vector u(t) at every time step t (t ∈ 0, T − 1) is zonotope, i.e. it can be represented as  P(t) = c2 , hl1 , . . . , lq i , q ∈ N. Thus, the following proposition holds. Proposition 1. Let initial set X0 and control constraint set P(t) of the linear discrete-time dynamical system (1) be zonotopes. Then the exact reachable set G 0, X0 , T of this dynamical system at final (terminal) time T is also zonotope. Proof. Since zonotopes are the special  class of convex polytopes, it also have a evolutionary property, that is G 0, X0 , t + 1  = G t, X(t), t + 1 and X(t) = G 0, X0 , t), t ∈ 1, T − 1. As it shown above, the exact reachable set G 0, X0 , T at final (terminal) time T can be determined by recurrent sequence of computation of transitional  reachable sets X(t + 1) = G t, X(t), t + 1 , t ∈ 0, T − 1. Thus, it is necessary to show that transitional reachable  sets G t, X(t), t + 1 , t ∈ 0, T − 1 are zonotopes. For t = 0 by definition of reachable set we have   X(1) = G 0, X0 , 1 = x(1) : x(1) = Ax(0) + Bu(0), x(0) ∈ X0 , u(0) ∈ P(0) . Taking into account Assumptions 1 and 2, we can consider this reachable set as Minkowski sum of two sets:  X(1) = G 0, X0 , 1 = AX0 + BP(0) =   = Ac1 , hAr1 , . . . , Arp i + Bc2 , hBl1 , . . . , Blq i =  = Ac1 + Bc2 , hAr1 , . . . , Arp , Bl1 , . . . , Blq i , p ∈ N, q ∈ N. It is clear that X(1) is a zonotope. Accordingly, using X(1) as the initial set it is easy to find that X(2) is also a zonotope. By induction, using zonotope X(T − 1) as initial set we define that the set G T − 1, X(T − 1), T is a zonotope. This completes the proof.  Note that to obtain all vertices of some zonotope Z we assume that the generators are of form ri = [0, bi ] for all i = 1, . . . , m. For each i we just select 0 or bi and compute the sum. Each vertex of Z is such a sum. There are 2m sums that constitute this highly redundant V-representation. Then, one can remove all redundant ones using linear programming methods (e.g. [Pa95], [Sh93]). This is a simple but very inefficient way to compute V-polytope from zonotope. It can be employed only for low dimensions of state space and little number of time steps. However, there is a much more efficient algorithm [AF96] which is polynomial in the size of both input and output. Now, we ready to present an algorithm of reachable set computation using zonotopes. 5 Algorithm of reachable set computation An algorithm presented in this section is a modification of one in ([Sh93]), and it consists in using of zonotopes instead of V-polytopes. The obtaining of the algorithm is based on Proposition 1 and the results of previous section. The purpose of the following recurrent algorithm is to compute reachable set of system (1)–(3) on one step forward. Algorithm 1 (Reachable set computation). Phase 1: Generating the set of states described by zonotope:  ZX (t) = c1 , hr1 , . . . , rn i , where c1 – center of zonotope X(t), hr1 , . . . , rn i – generators of X(t), t ∈ 0, T − 1. 107 Phase 2: Generating the set of feasible controls described by zonotope:  ZP (t) = c2 , hl1 , . . . , lp i , where c2 – center of zonotope P(t), hl1 , . . . , lp i – generators of P(t), t ∈ 0, T − 1. Phase 3: Linear transformation of zonotopes:  ẐX (t) = AZX (t) = Ac1 , hAr1 , . . . , Arn i ,  ẐP (t) = BZP (t) = Ac2 , hBl1 , . . . , Blp i . Phase 4: Determination of the Minkowski sum:  X(t + 1) = ẐX (t) + ẐP (t) = Ac1 + Ac2 , hAr1 , . . . , Arn , Bl1 , . . . , Blp i , which describes reachable set of system (1)–(3) at moment t + 1. It is necessary to note that if matrix A in (1) is not a unity matrix, then the number of reachable set generators will grow linearly. But when one need to compute V-representation of reachable set from zonotope representation, the linear programming problem of convex hull computation can expend too much time. Thus, computational capability of PC constrains in a manner the number of time steps T and dimensionality of a system. 6 Modeling example The obtained reachable set computation algorithm 1 has been realized in MATLAB R2014a. We are going to show the performance of proposed algorithm in comparison to exact reachable set computation algorithm using polytopes at the numerical simulation example of launch vehicle two-tank sixteen-step discrete-time dynamical system in R4 , presented in Section 2. Example. On the given time interval 0, 16 the dynamic of system (1)–(3), which reproduce the launch vehicle propulsion system work, is considered. The state vector is: x(t) = x1 (t), x2 (t), x3 (t), x4 (t) ∈ R4 , t ∈ 0, 16.  Control of the system is a scalar u(t) ∈ R, t ∈ 0, 15. Matrices A and B in (1) are equal to:     1 0 0 0 4 −6.5 1 0 0  0 A=  0 , B =  . 0 1 0 2 0 0 −6.5 1 0 The initial set of states X0 is a zonotope with:           1150 20 0 0 0 119600 0 500 0  0   430  , r1 =  0  , r2 =  0  , r3 = 10 , r4 =  0  . c1 =           44720 0 0 0 300 Since u(t) 6 1, the feasible control set P(t) is a zonotope with:     0 1 0 1 c2 = 0 , l1 = 1 .    0 1 Figure 3 shows projections of reachable sets on x1 –x2 , and x3 –x4 axes, correspondingly. The black dots on Fig- ure 3 illustrate the projection of reachable sets computed by algorithm using V-Polytopes. Also, we demonstrate projections of reachable sets in R3 computed using presented algorithm. These projections in R3 are shown in Figure 4. 108 Figure 3: The projections of reachable sets in R2 using algorithms for zonotopes and polytopes Figure 4: The projections of reachable sets in R3 109 The simulation was realized on an Intel Core i5 3.4 GHz with 8 Gb RAM. The time of simulation is less than 0.5 sec for both algorithm using zonotopes and algorithm using V-polytopes. It is necessary to note that for systems with dimensions higher than 8 and time steps higher than, say, 50, both algorithms becomes inefficient since the simplex-method is used. 7 Conclusion In this paper, we presented an algorithm for computation of exact reachable set of linear discrete-time systems. That algorithm is a modification of one presented in [Sh93]. The modification uses zonotopes, which allow an efficient implementation of the algorithm to fuel consumption terminal control problem, since time steps are always low and the sets of constraints may be considered as zonotopes. The use of zonotopes for computation of exact reachable sets on large time intervals is possible only using some approximation procedures (e.g. [Co03], [Gi05]), but there will be a loss in accuracy of reachable set. However, there are methods allow us to find all vertices of zonotope in polynomial time in the size of both input and output [AF96]. Our future work should focus on several points. First, using reachable set computation algorithm we should solve open-loop fuel consumption optimal terminal control problem of liquid-propellant vehicle. Second, it is necessary to extend the algorithm to handle linear discrete-time system with external uncertain disturbances. The third and the most important point is to solve closed-loop fuel consumption optimal terminal control problem of liquid-propellant vehicle. 7.0.1 Acknowledgements This paper is performed with the assistance of Russian Foundation for Basic Research (project №15-01-02368). References [AI13] A.Y. Andrienko, V.P. Ivanov. Theoretical and Practical Problems in On-Board Terminal Systems Design for Liquid-Propellant Carrier Rockets. Automation and Remote Control, 74(3):401–412, 2013. [ASB10] M. Althoff, O. Stursberg, M. Buss. Computing Reachable Sets of Hybrid Systems Using a Combination of Zonotopes and Polytopes. Nonlinear Analysis: Hybrid Systems, 4:233–249, 2010. [AF96] D. Avis, K. Fukuda. Reverse Search for Enumeration. Discrete Appl. Math., 65:21–46, 1996. [Co03] C. Combastel. A State Bounding Observer Based on Zonotopes. European Control Conference (ECC), 2589–2594, 2003. [Fu04] K. Fukuda. From the Zonotope Construction to the Minkowski Addition of Convex Polytopes. Journal of Symbolic Computation, 38:1261–1272, 2004. [Gi05] A. Girard. Reachability of Uncertain Linear Systems Using Zonotopes. Hybrid Systems: Computation and Control, 3414:291–305, 2005. [Pa95] P.M. Pardalos, Y. Li, W.W. Hager. Linear Programming Approaches to the Convex Hull Prob- lem in Rm . Computers & Mathematics with Applications, Volume 29, №7 : 23–29, 1995. [Sh93] V.A. Tyulyukin, A.F. Shorikov. Algorithm for Handling the Terminal Control Problem for Linear Discrete System. Automation and Remote Control, 4:115–127, 1993. [SK16-1] A.F. Shorikov, V.I. Kalev. Linear Discrete-Time Dynamical Model Forming for Solving Optimal Ter- minal Fuel Consumption Problem of Launch Vehicle. Informatsionnye tekhnologii i sistemy: sbornik trudov konferentsii, 61–66, 2016.(in Russian) [SK16-2] A.F. Shorikov, V.I. Kalev. Fuel Consumption Terminal Control Problem Statement for Liquid- Propellant Rockets. Izvestiya vuz. Fizika, 45–48, 2016. (in Russian) [Zi98] G.M. Ziegler. Lectures on Polytopes. New York, Springer Verlag, 1998. 110