Cost-limited reachability: a new problem in reachability analysis
CCost-limited reachability: a new problem inreachability analysis
Wei Liao ∗ a,b , Xiaohui Wei † a,b , and Jizhou Lai ca Key laboratory of Fundamental Science for National Defense-Advanced Design Technology of Flight Vehicle, NanjingUniversity of Aeronautics and Astronautics, Nanjing, Jiangsu, China b State Key Laboratory of Mechanics and Control of Mechanical Structures, Nanjing University of Aeronautics andAstronautics, Nanjing, Jiangsu, China c Navigation Research Center, Nanjing University of Aeronautics and Astronautics, Nanjing, Jiangsu, China
January 26, 2021
Abstract
In this paper, we generalize the definition of backward reachable tube. The classic definition ofbackward reachable tube is a set of system state that can be driven into the target set within a giventime horizon. Sometimes, the concern of researchers is not only the time consumption, but someother forms of cost of driving the system state toward the target set. Under this background, thedefinition of cost-limited backward reachable tube is put forward in this paper, where the cost is thetime integral of a running cost. And the running cost is a scalar function of system state and controlinput. A method to compute the cost-limited backward reachable tube is proposed. In this method,a cost-limited backward reachable tube is characterized by a non-zero level set of a value function,which is approximated using recursion and interpolation. At the end of this paper, some examplesare taken to illustrate the validity and accuracy of the proposed method.
Reachability analyses is an important method to study the behavior of nonlinear control systems. Inreachability problems, one specifies a target set in the state space and a time horizon and then seeksto find the backward reachable tube (BRT) of the target set [23]. However, computing the BRT is achallenging task. Due to the diversity of system states and control inputs, verifying each state one byone always consumes a lot of computation time [6, 4]. This makes the approaches based on simulationinsufficient, and thus formal verification methods are needed. Unfortunately, except for some solvableproblems [9], to exactly compute these sets is generally intractable. Therefore, the common practice is tocompute their approximate expressions. Numerous techniques have been proposed in the past decades,which can be divided into two categories: the Lagrangian methods [24] and the methods based on statespace discretization. Lagrangian methods mainly include ellipsoidal methods [18, 33, 19], polyhedralmethods [10, 1, 7] and support vector machines [24], which are able to solve high-dimensional problemsbut have strict requirements on the form of dynamic systems and are primarily used to solve linearproblems. The main advantage of the state space discretization-based approaches is its less stringentrequirements on the form of the dynamic system, and it can solve nonlinear problems. This category ∗ [email protected] † Corresponding author: wei [email protected] a r X i v : . [ ee ss . S Y ] J a n f methods mainly includes level set method [23, 6, 25] and distance fields over grids (DFOG) method[11, 2, 3], among which the level set method is more widely used in engineering.In level set method, the Hamilton-Jacobi (HJ) partial differential equation (PDE) about a valuefunction is numerically solved, and the RBS can then be expressed as the zero level set of solution of theHJ PDE [23]. Based on this principle, some mature toolboxes [27, 26, 14, 12] have been developed andhave greatly promoted the development of computational techniques for RBS and solved many practicalengineering problems, such as flight control systems [32, 29, 34], ground traffic management systems[8, 20], air traffic management systems [5], etc.Although the level set method has developed into a mature, reliable and commonly used tool overthe years, there are still some shortcomings in this method, the most important of which may be thedifficulty of handling the dissipation function [27]. The dissipation function is an optimization problemduring the process of solving the viscosity solution of the HJ PDE. To our best knowledge, although thisoptimization problem can be reduced in affine nonlinear systems, there is no efficient method to solve it inmore general systems. Literature [27] suggests that the instability can be avoided by overapproximatingthe dissipation function. However, this may lead to relatively large calculation errors. Furthermore, themechanism of the level set method determines the necessity to store the solutions of the HJ PDE atdifferent moments to store the BRTs within different time horizons.To some extent, a reachability problem can be regarded as an extension of an optimal control problem.The classical definition of reachable tube focuses on whether a state can be driven into the target state ina limited time, it has a deep logical interrelationship with time optimal control problems. However, timeoptimal control problem is merely a special form of optimal control problem. In a general form of optimalcontrol problem, the performance index always consists of a Lagrangian [31], which is the time integralof a running cost. Therefore, it is quite necessary to develope a novel framework to discuss the ability ofreaching the target set within an admissible performance index. At the same time, such a framework alsohas a great engineering significance. For example, fuel consumption is a common performance index insome engineering practices, researchers may concern whether the target set can be reached before runningout of fuel.To overcome the shortcoming of level set method and discuss the reachability under different formsof running cost, this paper proposes a more generalized definition of BRT and a novel framework to solvereachability problems. In summary, the main contributions of this paper are as follows:(1) We put forward the definition of cost-limited backward reachable tube (CBRT). In our research,a running cost function can be specified for the system, the CBRT is a set of states that can bedriven into the target set before the performance index increasing to the given admissible cost.(2) We propose a method to compute the CBRT (and the BRT). In this method, the CBRT can bedescribed as a level set (non-zero level set) of a value function. This method differs from level setmethod in that the value function is obtained not by solving the HJ PDE but by recursion andinterpolation. Such a mechanism avoids solving PDE and thus can effectively avoid the computa-tion of the dissipation function. In addition, the CBRTs under different admissible costs can berepresented as different level sets of the same value function. This allows for significant storagespace saving.The structure of this paper is as follows. Section II introduces the definitions of BRT and CBRT.Section III describes the proposed method in detail. Some numerical examples are given in Section IVto illustrate the validity of the proposed method. The results are summarized in Section V.2 Reachability and optimal control
Consider a continuous time control system with fully observable state and non-deterministic boundeddisturbance: ˙ s = f ( s, a, d ) (1)where s ∈ R n is the system state, a ∈ A ⊂ R n and d ∈ D ⊂ R n are refered to as the controller and thedisturbance. The function f : R n × A × D is Lipschitz continuous and bounded. Let A and D denotethe set of lebesgue measurable functions from the interval [0 , ∞ ) to A and D respectively. Then, giventhe state s at time t , a ( . ) ∈ A , d ( . ) ∈ D , the evolution of system (1) in time interval [ t , t ] can bedenoted as a continuous trajectory φ t t ( ., s , a, d ) : [ t , t ] → R n and φ t t ( t , s , a, d ) = s . Given a targetset K and a time horizon ¯ T , the definition of backward reachable tube (BRT) is [23, 25]: Definition 1.
Backward reachable tube (BRT): R K ( ¯ T ) = (cid:8) s ∈ R n |∃ t ∈ [0 , ¯ T ] , ∀ d ( . ) ∈ D , ∃ a ( . ) ∈ A ,φ ¯ T ( t, s , a, d ) ∈ K (cid:9) (2)In some literatures [23], the reachability problem is expressed as the following optimal control problem: V ( s , t ) = inf a ( . ) ∈ A sup d ( . ) ∈ D min τ ∈ [ t, ¯ T ] l (cid:16) φ ¯ T ( τ, s , a, d ) (cid:17) s.t. V ( s , ¯ T ) = l ( s ) (3)where l ( . ) is bounded and Lipschitz continuous, and satisfies K = { s ∈ R n | l ( s ) ≤ } . The BRT canthen be characterized by the zero level set of V ( ., R K ( ¯ T ) = { s ∈ R n | V ( s , ≤ } (4)In this section, we derive the representation of the BRT from another form of optimal control problem,and generalize it to CBRT.Consider a case where the control input a ( . ) is chosen to enable the system state to reach the targetset in as short a time as possible, and the disturbance d ( . ) is chosen to avoid or delay as much as possiblethe entry of the system state into the target set. If a trajectory can enter the target set within timehorizon ¯ T in this case, then the initial state of this trajectory must belong to the BRT. Thus, a value3unction W ( . ) : R n → R can be constructed as Equ. (5). W ( s ) = min a ( . ) max d ( . ) t f s.t. ˙ s ( t ) = f ( s ( t ) , a ( t ) , d ( t )) ∀ t ∈ [0 , t f ] s (0) = s a ( t ) ∈ A ∀ t ∈ [0 , t f ] d ( t ) ∈ D ∀ t ∈ [0 , t f ] s ( t f ) ∈ K , if ∀ d ( . ) ∈ D ∃ a ( . ) ∈ A ∃ t ∈ [0 , ∞ ) φ ∞ ( t, s , a, d ) ∈ K ∞ , otherwise (5)and the BRT can be characterized by the level set of ¯ T of W ( . ), i.e.: R K ( ¯ T ) = (cid:8) s ∈ R n | W ( s ) ≤ ¯ T (cid:9) (6) As described earlier, for a general optimal control problem, the performance index always consists of aLagrangian. Denote the running cost as c ( ., . ) : R n × A → R , and assume that: Assumption 1. min s ∈ R ,a ∈A c ( s, a ) = γ holds, where γ is a positive real number. The performance index of the evolution initialized from s at time t under control input a ( . ) anddisturbance d ( . ) in time interval [ t , t ] is denoted as: J t t ( s , a, d ) = (cid:90) t t c (cid:0) φ t t ( t, s , a, d ) , a ( t ) (cid:1) dt (7)Then given a target set K and an admissible cost ¯ J , the CBRT can be defined: Definition 2.
Cost-limited backward reachable tube (CBRT): R cK ( ¯ J ) = (cid:8) s ∈ R n |∃ t ∈ [0 , ∞ ) , ∀ d ( . ) ∈ D , ∃ a ( . ) ∈ A ,φ ∞ ( t, s , a, d ) ∈ K ∧ J t ( s , a, d ) ≤ ¯ J (cid:9) (8)Informally, under Assumption 1, the performance index J t ( s , a, d ) increases with the increasing of t , R cK ( ¯ J ) is a set of states that can be transfered into the target set K under any disturbance d ( . ) beforethe performance index increasing to the given admissible cost ¯ J .Consider a case where the controller aims to transfer the system state to the target set with the leastpossible cost, and the disturbance aims to avoid the entry of the system state into the target set or tomaximize the cost during state transition. If a trajectory can enter the target set before the performanceindex increasing to the given admissible cost ¯ J , then the initial state of this trajectory must belong to4 cK ( ¯ J ). Similar to Equ. (5), a value function can be constructed as follows: W c ( s ) = min a ( . ) max d ( . ) (cid:90) t f c ( s ( t ) , a ( t )) dt s.t. ˙ s ( t ) = f ( s ( t ) , a ( t ) , d ( t )) ∀ t ∈ [0 , t f ] s (0) = s a ( t ) ∈ A ∀ t ∈ [0 , t f ] d ( t ) ∈ D ∀ t ∈ [0 , t f ] s ( t f ) ∈ K , if ∀ d ( . ) ∈ D ∃ a ( . ) ∈ A ∃ t ∈ [0 , ∞ ) φ ∞ ( t, s , a, d ) ∈ K ∞ , otherwise (9)and the CBRT can be characterized by the level set of ¯ J of W c ( . ), i.e.: R cK ( ¯ J ) = (cid:8) s ∈ R n | W c ( s ) ≤ ¯ J (cid:9) (10) Remark 1. If c ( s, a ) ≡ , then (cid:90) t f c ( s ( t ) , a ( t )) dt = t f , the general optimal control problem in Equ. (9)degenerate into the time optimal control problem in Equ. (5). Similarly, the CBRT R cK ( ¯ J ) is degeneratedinto the BRT R K ( ¯ J ) . Generally speaking, the relationship between optimal control and reachability can be intuitively ex-pressed with Fig. 1.
Optimal controlproblemsCBRS Time optimal controlproblemsBRS
Figure 1: The relationship between optimal control and reachability.
The optimal control problem in Equ. (9) is an optimal control problem with terminal state constraint.In order to convert it into an optimal control problem without terminal state constraint to facilitate theconstruction of the recursive formula of the value function, we shall need to define a modified dynamicsystem: ˙ s = ˆ f ( s, a, d ) = f ( s, a, d ) , s / ∈ K , s ∈ K (11)and a modified running cost function: ˆ c ( s, a ) = c ( s, a ) , s / ∈ K , s ∈ K (12)5iven the state s at time t , a ( . ) ∈ A , d ( . ) ∈ D , the evolution of the modified system (11) in timeinterval [ t , t ] can also be expressed as a continuous trajectory ˆ φ t t ( ., s , a, d ) : [ t , t ] → R n , and themodified performance index is:ˆ J t t ( s , a, d ) = (cid:90) t t ˆ c (cid:16) ˆ φ t t ( t, s , a, d ) , a ( t ) (cid:17) dt (13) Remark 2.
As long as the trajectory ˆ φ t t ( ., s , a, d ) evolves outside the target set K , it is the same astrajectory φ t t ( ., s , a, d ) , and the modified running cost function ˆ c ( s, a ) is equal to the c ( s, a ) . When thetrajectory ˆ φ t t ( ., s , a, d ) touches the border of K , then it stays on the border under the dynamics (11),and the modified running cost function ˆ c ( s, a ) is equal to zero. Then, a modified value function can be defined based on an optimal control problem without terminalstate constraint: (cid:99) W ct f ( s ) = min a ( . ) max d ( . ) (cid:90) t f ˆ c ( s ( t ) , a ( t )) dt s.t. ˙ s ( t ) = ˆ f ( s ( t ) , a ( t ) , d ( t )) ∀ t ∈ [0 , t f ] s (0) = s a ( t ) ∈ A ∀ t ∈ [0 , t f ] d ( t ) ∈ D ∀ t ∈ [0 , t f ] (14)With expression (13) introduced, the above equation can be simplified to the following form: (cid:99) W c ( s ) = min a ( . ) ∈ A max d ( . ) ∈ D ˆ J t f ( s , a, d ) (15)Then we may deduce the following result: Theorem 1.
For any ˜ J < γt f and any s ∈ { s ∈ R n | (cid:99) W ct f ( s ) ≤ ˜ J } , (cid:99) W ct f ( s ) = W c ( s ) holds. Proof:
According
Remark 2 , If the trajectory ˆ φ t f ( ., s , a, d ) cannot reach the target state K withintime horizon t f , then its performance index satisfies:ˆ J t f ( s , a, d ) ≥ (cid:90) t f min s ∈ R ,a ∈A c ( s, a ) dt = γt f (16)Therefore, (cid:99) W ct f ( s ) < ˜ J ≤ γt f implies that for any d ( . ) ∈ D , there exists a a ( . ) ∈ A and a τ ∈ [0 , t f )such that ˆ φ t f ( τ, s , a, d ) ∈ K . In such case, the following equations hold: (cid:99) W ct f ( s ) = min a ( . ) ∈ A max d ( . ) ∈ D (cid:90) t f ˆ c (cid:16) ˆ φ t f ( t, s , a, d ) , a ( t ) (cid:17) dt = min a ( . ) ∈ A max d ( . ) ∈ D (cid:90) τ ˆ c (cid:16) ˆ φ t f ( t, s , a, d ) , a ( t ) (cid:17) dt + (cid:90) t f τ dt = min a ( . ) ∈ A max d ( . ) ∈ D (cid:90) τ c (cid:16) φ t f ( t, s , a, d ) , a ( t ) (cid:17) dt s.t. ˆ φ t f ( τ, s , a, d ) ∈ K = W c ( s ) (17) (cid:3) Theorem 1 indicates that, given a modified value function (cid:99) W ct f ( . ), for any ¯ J < γt f , the CBRT R cK ( ¯ J )6an be represented by: R cK ( ¯ J ) = (cid:110) s ∈ R n | (cid:99) W ct f ( s ) ≤ ¯ J (cid:111) (18) According to Equ. (18), the crucial point in characterizing the CBRT is to compute the modified valuefunctions. For any τ ∈ [0 , t f ], the following equations hold: (cid:99) W ct f ( s ) = min a ( . ) ∈ A max d ( . ) ∈ D ˆ J t f ( s , a, d )= min a ( . ) ∈ A max d ( . ) ∈ D (cid:20) ˆ J τ ( s , a, d )+ min a (cid:48) ( . ) ∈ A max d (cid:48) ( . ) ∈ D ˆ J t f τ ( s (cid:48) , a (cid:48) , d (cid:48) ) (cid:21) = min a ( . ) ∈ A max d ( . ) ∈ D (cid:20) ˆ J τ ( s , a, d ) + (cid:99) W ct f − τ ( s (cid:48) ) (cid:21) (19)Divide the time interval into parts with length ∆ t , if ∆ t is small enough, a ( . ) and d ( . ) can be regardedas a constant in interval [ k ∆ t, ( k + 1)∆ t ] for any k ∈ N . And the modified system in Equ. (11) can beconverted into the discretized form: s (( k + 1)∆ t ) = (cid:98) F ( s ( k ∆ t ) , a ( k ∆ t ) , d ( k ∆ t )) (20)Note that, A simple approach to derive Equ. (20) is Euler’s method. For higher accuracy, a high orderRunge-Kutta method can also be used.According to trapezoid rule, the integral of the modified running cost function over the interval[ k ∆ t, ( k + 1)∆ t ] is approximated using the following equation: (cid:98) C ( s ( k ∆ t ) , a ) = ˆ c ( s ( k ∆ t ) , a ) + ˆ c ( F ( s ( k ∆ t ) , a, d ) , a )2 ∆ t (21)The recursive formula of the modified value function can then be obtained: (cid:99) W c ( s ) = 0 (cid:99) W c ( k +1)∆ t ( s ) = min a ∈A max d ∈D (cid:20) (cid:98) C ( s , a )+ (cid:99) W ck ∆ t (cid:16) (cid:98) F ( s , a, d ) (cid:17) (cid:21) (22)Given an admissible cost ¯ J , the required number of recursions can be determined by the followingtheorem. Theorem 2.
Given an admissible cost ¯ J , it takes no more than h = Int (cid:16) ¯ Jγ ∆ t (cid:17) + 1 times of recursionsto compute R cK ( ¯ J ) , where Int( x ) gives the integer closest to x . Proof:
The value function (cid:99) W ch ∆ t ( . ) can be obtained after h times of recursions. Since h ∆ tγ > ¯ J ,according to Theorem 1 , for any s ∈ { s ∈ R n | (cid:99) W ch ∆ t ( s ) ≤ ¯ J } , (cid:99) W ch ∆ t ( s ) = W c ( s ) holds. Therefore, R cK ( ¯ J ) = (cid:8) s ∈ R n | W c ( s ) ≤ ¯ J (cid:9) = (cid:110) s ∈ R n | (cid:99) W ch ∆ t ( s ) ≤ ¯ J (cid:111) (cid:3) The proposed method involves the solution of the min-max problem, which is in general hard to solve.In the current research, we use the sets of discrete points uniformly distributed in A and D when solvingthe min-max problem. Although a sufficient number of discrete points are required to accurately estimate7he value of min-max problem, in practice we find that in most examples the number of points does notsignificantly affect the computational results. In some special cases, such as the computation of the BRTsof affine nonlinear systems, only the endpoints of the intervals in A and D need to be checked. Normally, for any k ∈ Z + , it is difficult to compute the analytic form of (cid:99) W ck ∆ t ( . ). In the current research,a rectangular domain in the state space is specified as the computational domain and is discretized into aCartesian grid structure. The value of (cid:99) W ck ∆ t ( . ) at each grid point is stored in an array whose number ofdimensions is the same as that of the dynamic system. Take a three-dimensional system as an example,where s = [ x, y, z ] T . In this case, the function (cid:99) W ck ∆ t ( . ) can be represented by a trilinear interpolation[30] of the aforementioned array. The pseudocode of the proposed method is shown in Algorithm 1. Algorithm 1
Computation of modified value functions Input:
Recursion number h , number of grid points N x , N y , N z , computational domain Ω = [ x d , x u ] × [ y d , y u ] × [ z d , z u ]; Construct the modified system in discretized form (20); Construct the modified running cost in discretized form (21); Let W and W (cid:48) be two N x × N y × N z arrays and initialize them to ; δx ← x u − x d N x − , δy ← y u − y d N y − , δz ← z u − z d N z − for k ← , ..., h do Construct the trilinear interpolation function (cid:99) W ( . ) using W ; for i x ← , ..., N x − do for i y ← , ..., N y − do for i z ← , ..., N z − do s ← x d + i x δxy d + i y δyz d + i z δz W (cid:48) [ i x ][ i y ][ i z ] ← min a ∈A max d ∈D (cid:20) (cid:98) C ( s , a )+ (cid:99) W ck ∆ t (cid:16) (cid:98) F ( s , a, d ) (cid:17) (cid:21) end for end for end for Copy W (cid:48) to W ; end for Return (cid:99) W ( . ); In this section, two examples are taken to illustrate the validity of the proposed method. The firstexample shows the BRT of a two-dimensional linear system, which can be analytically solved, to analysisthe effects of parameters in
Algorithm 1 on computational accuracy. The second example, on thelongitudinal flight envelope of an aircraft, is a practical application of CBRT.8 .1 Two-dimensional linear system example
Consider the following system: ˙ s = (cid:34) ˙ x ˙ y (cid:35) = (cid:34) d − x + a (cid:35) (23)where s = [ x, y ] T ∈ R is the system state, a ∈ A = [ − ,
1] is the control input, d ∈ D ∈ [ − ,
1] is thedisturbance. The target set K = (cid:8) [ x, y ] T || y | ≤ . (cid:9) , the time horizon ¯ T = 1. To compute the BRT, therunning cost function is set as c ( s, a ) ≡
1. The analytical solution of this problem is: R ∗ K ( ¯ T ) = K ∪ { [ x, y ] T | ≤ y ≤ x + 1 , x ≥ }∪{ [ x, y ] T | ≥ y ≥ x − , x ≤ }∪{ [ x, y ] T | ≤ y ≤ x + x + 1 , − ≤ x ≤ }∪{ [ x, y ] T | ≥ y ≥ − x + x − , ≤ x ≤ } (24)In order to analysis the convergence of the proposed method, we compute the BRT with several sets ofparameters. Fig. 2 shows the results of our method. The number of recursions in each of these examplesis h = Int (cid:0) t (cid:1) +1. Since system (23) is a two-dimensional system, the bilinear interpolation [ ? ] is appliedin our method.Figure 2: Computational results under different parameters (the analytical solution is expressed in gray).9o quantify the computational errors under different parameters, the relative volume errors is quanti-fied by Jaccard index [13] between R ∗ K ( ¯ T ) and R K ( ¯ T ) = { s ∈ R n | (cid:99) W ch ∆ t ( s ) ≤ } , i.e, 1 − | R ∗ K ( ¯ T ) ∩ R K ( ¯ T ) || R ∗ K ( ¯ T ) ∪ R K ( ¯ T ) | .Table 1 summarizes these relative volume errors.Table 1: Relative error∆ t G r i dnu m b e r
51 0.0305 0.0304 0.0304 0.0304101 0.0243 0.0223 0.0207 0.0192151 0.0158 0.0155 0.0157 0.0167201 0.0211 0.0188 0.0188 0.0178Counter-intuitively, it is not the case that the smaller the time step, the more accurate the result.This may be due to the accumulated interpolation error, which decreases with the decreasing number ofrequired recursions.Moreover, in our method, the BRTs under different time horizons ¯ T (or the CBRTs under differentadmissible costs ¯ J ) can be expressed as different level sets of the same value function, see Fig. 3.Figure 3: The BRTs under different time horizons. Longitudinal motion is the subsystem in the aircraft’s plane of symmetry ( x − z plane) and is characterizedby the following states: velocity v , attack angle α , pitch rate q , and pitch angle θ . These variables aredescribed in the ground coordinate system, see Fig. 4. The dynamic equation of the longitudinal motioncan be written as Equ. (25) [29]:˙ s = (cid:104) ˙ v, ˙ α, ˙ q, ˙ θ (cid:105) T = m (cid:0) P cos α − ρv SC D ( α, δ e ) − mg sin( θ − α ) (cid:1) q + mv (cid:0) − P sin α − ρv SC L ( α, δ e ) + mg cos( θ − α ) (cid:1) I y (cid:0) ρv S ¯ bC m ( α, q, δ e ) + d (cid:1) q (25)where s = [ v, α, q, θ ] T is system state, δ e ∈ [ − . , .
3] is the elevator deflection and is regarded as thecontrol input, d ∈ [ − . × , . × ] is disturbance moment. and P is the thrust, m is the mass of10igure 4: Longitudinal motion of an aircraft.the aircraft, ρ is the air density, I y is the inertia moment around y axis, C L , C D , C m are the aerodynamicderivatives. In this example, it is assumed that the thrust P is varied to keep the velocity v be a constant,that is: ˙ v = 0 = ⇒ P = 1cos α (cid:20) ρv SC D ( α, δ e ) + mg sin( θ − α ) (cid:21) (26)Under this assumption, the state space degenerates to a three-dimensional space, i.e., s = [ α, q, θ ] T . Theaerodynamic model is constructed as the following equations: C L = C L + C Lα α + C Lδ e δ e C D = C D + C Dα α + C Dα α + C Dδ e δ e C m = C m + C mα α + C mq q ¯ b v + C mδ e δ e (27)The parameters in Equ. (25) and Equ. (27) are summarized in Table 2. In this example, the task is toTable 2: Physical parameters of the reference aircraft Parameter Value Parameter Value m I y S
524 ¯ b ρ C L C Lα v C Lδe C D C Dα C Dα C mδe -1.2 C Dδe C m C mα -0.2 C mq -1drive the state trajectory into the target set K = [ − . , . × [ − . , . × [ − . , .
1] within the givenadmissible cost ¯ J = 1. The running cost is a weighted sum of time consumption and maneuver overloadfactor, i.e.: c ( s, δ e ) = 1 + λ G (cid:107) G (cid:107) = 1 + λ G mg (cid:112) | F z | + | F x | (28)11n the preceding equation, λ G is the weight of the maneuver overload, F x and F z are the components ofthe summation of aerodynamic force and thrust in the horizontal and vertical direction respectively, i.e.: F x = P cos θ − ρv S [ C L sin( θ − α ) − C D cos( θ − α )] F z = P sin θ + 12 ρv S [ C L cos( θ − α ) − C D sin( θ − α )] (29)We consider two cases: λ G = 0 and λ G = 0 .
25. In the former case, the running cost c ( s, δ e ) ≡
1, suchthat the problem degenerates to the computation of BRT, which can be solved by level set method. Thesolver settings of our method for the BRT computation are: • Computational domain: Ω = [ − . , . × [ − . , . × [ − . , . • Number of grid points: 101 × × • Time step size: ∆ t = 0 . • Number of recursions: 101.The comparison between the results of the proposed method and level set method is shown in Fig. 5(The setups of the computational domain and the grid points’ number in level set method are the sameas those in our method). It can be seen that, the envelope of the BRT computed by the proposed methodalmost coincides with that computed by level set method.Figure 5: Computation results of the BRT (the result of our method is expressed in blue and the resultof the level set method is expressed in red). 12he later case is about the computation of CBRT. The solver settings of our method are the same asthose in the previous case. The computation results of CBRTs under different ¯ J are shown in Fig. 6.Figure 6: CBRTs under different ¯ J . An extension of backward reachable tube (BRT) refered to as cost-limited backward reachable tube(CBRT) is introduced in this paper. A CBRT is a set of state that can be driven into the target setbefore the performance index growing to a given admissible cost, and the performance is the time integralof a running cost function, which is a scalar function of system state and control input. In addition, amethod to compute the CBRT is proposed. In this method, a CBRT is characterized by a non-zero levelset of a value function which are approximated by recursion and interpolation.The primary weakness of the proposed method to compute CBRT, and many other methods to com-pute BRT, is the exponential growth of memory and computational cost as the system dimension increases[25, 17]. Some approaches have been proposed to mitigate these costs, such as projecting the BRT ofa high-dimensional system into a collection of lower dimensional subspaces [6, 28, 21] and exploitingthe structure of the system using the principle of timescale separation and solving these problems in asequential manner [15, 16, 22].Another open issue is mechanism for determining the time step size. As in the first example, a smallertime step does not necessarily lead to a more accurate computational result. These will be considered in13ur future work.
References [1] Eugene Asarin, Olivier Bournez, and Thao Dang. Approximate reachability analysis of piecewise-linear dynamical systems.
Approximate Reachability Analysis of Piecewise-Linear Dynamical Sys-tems , 1790/2000:21–31, 05 2000.[2] Robert Baier, Christof Buskens, I. Chahma, and Matthias Gerdts. Approximation of reachablesetsby direct solution methods for optimal control problems.
Optimization Methods & Software - OPTIMMETHOD SOFTW , 22:433–452, 06 2007.[3] Robert Baier, Matthias Gerdts, and Ilaria Xausa. Approximation of reachable sets using optimalcontrol algorithms.
Numerical Algebra, Control and Optimization , 3, 09 2013.[4] S. Bansal, M. Chen, S. Herbert, and C. J. Tomlin. Hamilton-jacobi reachability: A brief overviewand recent advances. In , pages2242–2253, 2017.[5] Somil Bansal, Mo Chen, Ken Tanabe, and Claire J Tomlin. Provably safe and scalable multivehicletrajectory planning.
IEEE Transactions on Control Systems Technology , 2020.[6] M. Chen, S. L. Herbert, M. S. Vashishtha, S. Bansal, and C. J. Tomlin. Decomposition of reach-able sets and tubes for a class of nonlinear systems.
IEEE Transactions on Automatic Control ,63(11):3675–3688, 2018.[7] Alongkrit Chutinan and Bruce H Krogh. Computational techniques for hybrid system verification.
IEEE transactions on automatic control , 48(1):64–75, 2003.[8] Aparna Dhinakaran, Mo Chen, Glen Chou, Jennifer C Shih, and Claire J Tomlin. A hybrid frameworkfor multi-vehicle collision avoidance. In , pages 2979–2984. IEEE, 2017.[9] T. Gan, M. Chen, Y. Li, B. Xia, and N. Zhan. Reachability analysis for solvable dynamical systems.
IEEE Transactions on Automatic Control , 63(7):2003–2018, 2018.[10] Mark R. Greenstreet. Verifying safety properties of differential equations. In Rajeev Alur andThomas A. Henzinger, editors,
Computer Aided Verification , pages 277–287, Berlin, Heidelberg,1996. Springer Berlin Heidelberg.[11] Roel Helsen, Erik-Jan van Kampen, Cornelis C. de Visser, and Qiping Chu. Distance-fields-over-grids method for aircraft envelope determination.
Journal of Guidance, Control, and Dynamics ,39(7):1470–1480, 2016.[12] Sylvia Herbert. helperoc. https://github.com/HJReachability/helperOC . Accessed 15 December2020.[13] Paul Jaccard. The distribution of the flora in the alpine zone. 1.
New phytologist , 11(2):37–50, 1912.[14] Tanabe K and Chen M. Beacls: Berkeley efficient api in c++ for level set methods. https://github.com/HJReachability/beacls . Accessed 15 December 2020.1415] I Kitsios and J Lygeros. Aerodynamic envelope computation for safe landing of the hl-20 personnellaunch vehicle using hybrid control. In
Proceedings of the 2005 IEEE International Symposium on,Mediterrean Conference on Control and Automation Intelligent Control, 2005. , pages 231–236. IEEE,2005.[16] I Kitsios and J Lygeros. Final glide-back envelope computation for reusable launch vehicle usingreachability. In
Proceedings of the 44th IEEE Conference on Decision and Control , pages 4059–4064.IEEE, 2005.[17] Ioannis Kitsios and John Lygeros.
Launch-Pad Abort Flight Envelope Computation for a PersonnelLaunch Vehicle Using Reachability .[18] Alexander B. Kurzhanski and Pravin Varaiya. Ellipsoidal techniques for reachability analysis. InNancy Lynch and Bruce H. Krogh, editors,
Hybrid Systems: Computation and Control , pages 202–214, Berlin, Heidelberg, 2000. Springer Berlin Heidelberg.[19] Alexander B Kurzhanski and Pravin Varaiya. Ellipsoidal techniques for reachability analysis. In
International Workshop on Hybrid Systems: Computation and Control , pages 202–214. Springer,2000.[20] Karen Leung, Edward Schmerling, Mengxuan Zhang, Mo Chen, John Talbot, J Christian Gerdes,and Marco Pavone. On infusing reachability-based safety assurance within planning frameworks forhuman–robot vehicle interactions.
The International Journal of Robotics Research , 39(10-11):1326–1345, 2020.[21] Anjian Li and Mo Chen. Guaranteed-safe approximate reachability via state dependency-baseddecomposition. In , pages 974–980. IEEE, 2020.[22] Thomas Lombaerts, Stefan Schuet, Kevin Wheeler, Diana M Acosta, and John Kaneshige. Safemaneuvering envelope estimation based on a physical approach. In
Aiaa guidance, navigation, andcontrol (gnc) conference , page 4618, 2013.[23] John Lygeros. On reachability and minimum cost optimal control.
Automatica , 40(6):917 – 927,2004.[24] John N Maidens, Shahab Kaynama, Ian M Mitchell, Meeko MK Oishi, and Guy A Dumont. La-grangian methods for approximating the viability kernel in high-dimensional systems.
Automatica ,49(7):2017–2029, 2013.[25] I. M. Mitchell, A. M. Bayen, and C. J. Tomlin. A time-dependent hamilton-jacobi formulation ofreachable sets for continuous dynamic games.
IEEE Transactions on Automatic Control , 50(7):947–957, 2005.[26] Ian Mitchell. A toolbox of level set methods. .Accessed 15 December 2020.[27] Ian M Mitchell. A toolbox of level set methods.
UBC Department of Computer Science TechnicalReport TR-2007-11 , 2007.[28] Ian M Mitchell. Scalable calculation of reach sets and tubes for nonlinear systems with terminalintegrators: a mixed implicit explicit formulation. In
Proceedings of the 14th international conferenceon Hybrid systems: computation and control , pages 103–112, 2011.1529] H. N. Nabi, T. Lombaerts, Y. Zhang, E. van Kampen, Q. P. Chu, and C. C. de Visser. Effects ofstructural failure on the safe flight envelope of aircraft.
Journal of Guidance, Control, and Dynamics ,41(6):1257–1275, 2018.[30] NASA. Trilinear interpolation. . Accessed December 4, 2020.[31] RWH Sargent. Optimal control.
Journal of Computational and Applied Mathematics , 124(1-2):361–371, 2000.[32] Jork Stapel, Coen De Visser, Erik-Jan Van Kampen, and Q. Chu. Efficient methods for flightenvelope estimation through reachability analysis. 01 2016.[33] Z. Xu, H. Su, P. Shi, R. Lu, and Z. Wu. Reachable set estimation for markovian jump neuralnetworks with time-varying delays.
IEEE Transactions on Cybernetics , 47(10):3208–3217, 2017.[34] Y. Zhang, C. C. de Visser, and Q. P. Chu. Database building and interpolation for an online safeflight envelope prediction system.