- Email: [email protected]
Chain of Set Inversion Problems Application to reachability analysis
Proceedings of the 20th World The International Federation of Congress Automatic Control Proceedings of the 20th World The International Federation of Congress Automatic Control Toulouse, France, July 9-14, 2017 Proceedings of the 20th World Congress The International Federation of Automatic Control Toulouse, France, July 9-14, 2017 Available online at www.sciencedirect.com The International of Automatic Control Toulouse, France,Federation July 9-14, 2017 Toulouse, France, July 9-14, 2017
ScienceDirect
IFAC PapersOnLine 50-1 Inversion (2017) 9266–9271 Problems Chain of Set Chain of Set Inversion Problems Chain of Set Inversion Problems Application to reachability analysis Chain of Set Inversion Problems Application to reachability Application to reachability analysis analysis Application to reachability analysis ∗ ∗∗ Benoit Desrochers Luc Jaulin ∗ ∗∗
Benoit Desrochers ∗ Luc Jaulin ∗∗ Benoit Desrochers ∗ Luc Jaulin ∗∗ Benoit DesrochersTechniques Luc Jaulin Direction G´en´erale de l’Armement, Navales Brest, BCRM Direction G´en´erale de l’Armement, Techniques Navales Brest, BCRM Brest, 29240 BrestTechniques Cedex 9, (e-mail: Direction G´en´ede rale de l’Armement, Navales Brest, BCRM Brest, 29240 BrestTechniques Cedex 9, (e-mail: Direction G´en´ede rale de l’Armement, Navales Brest, BCRM [email protected]) de Brest, 29240 Brest Cedex 9, (e-mail: [email protected]) ∗∗ de Brest, 29240 Brest Cedex 9, (e-mail: ENSTA-Bretagne, Lab-STICC, 2 rue Fran¸ c ois Verny, 29806 Brest, [email protected]) ∗∗ Lab-STICC, 2 rue Fran¸cois Verny, 29806 Brest, ∗∗ ENSTA-Bretagne, [email protected]) (e-mail: [email protected] )Verny, 29806 Brest, ENSTA-Bretagne, Lab-STICC, 2 rue Fran¸ c ois (e-mail: [email protected] ∗∗ ENSTA-Bretagne, Lab-STICC, 2 rue Fran¸cois))Verny, 29806 Brest, (e-mail: [email protected] (e-mail: [email protected] ) Abstract: This paper deals with the set inversion problem X = f −1 (Y) in the case where Abstract: This paper deals with the set inversion qproblem X = f −1 (Y) in the case where m −1 to f : Rnn → RThis on a parameter vector p ∈ Rqproblem which isX known a box [p]. (Y)beininside the case where Abstract: paper deals with the set inversion = f m depends f : Rn → RThis depends on a parameter vector p ∈ Rqproblem which isX known −1 to be inside a box [p]. m (Y) ininside the case where Abstract: paper deals with the set inversion = f We show that for a large class of problems, wep can obtain anis accurate approximation of the fWe : R → R depends on a parameter vector ∈ R which known to be a box show that for a large of problems, can obtain anis accurate of [p]. the n q f : R → Rm depends on class a parameter vectorwe ∈R which known toapproximation be a box [p]. solution set, without bisecting inofthe p-space. Top do this, symbolic methods areinside required to cast We show that for a large class problems, we can obtain an accurate approximation of the solution set, without bisecting in the p-space. To do this, symbolic methods are required to cast We show that for a into large can obtain an accurate approximation of cast the our initial problem a class chain ofproblems, set-inversion problems, the links of which have some nice solution set, without bisecting inofthe p-space. we To do this, symbolic methods are required to our initial problem into a chain of set-inversion problems, the links of which have some nice solution set, without bisecting in an the p-space. To we do this, symbolic methods required to cast properties with respect toa p. As application, consider thelinks problem ofare computing the set our initial problem into chain of set-inversion problems, the of which have some nice properties with respect toa p. As an application, we considerthe thelinks problem of computing thenice set our problem chain of set-inversion problems, have of allinitial initialwith states ofinto an uncertain discrete-time state system that reachofa which target set Y some in athe given properties respect to p. As an application, we consider the problem of computing set of all initial states of an uncertain discrete-time state system that reach a target set Y in a given properties with respect to p. As an application,state we consider the problem of computing the set time. of all initial states of an uncertain discrete-time system that reach a target set Y in a given time. of all initial states of an uncertain discrete-time state system that reach a target set Y in a given time. © 2017, IFAC (Internationalanalysis, Federation of Automatic Control) Hostingsystems; by Elsevier Ltd. All error rights reserved. time. Keywords: Reachability Interval analysis, Quantized Bounded Keywords: Reachability analysis, Interval analysis, Quantized systems; Bounded error identification, Keywords: Reachability analysis, Interval analysis, Quantized systems; Bounded error identification, Keywords: Reachability analysis, Interval analysis, Quantized systems; Bounded error identification, identification, 1. INTRODUCTION interval analysis are more and more used in the context of 1. INTRODUCTION interval analysis are more and more used in the context of nonlinearanalysis estimation Jauberthie et al. (2011) Raissi et al. 1. INTRODUCTION interval are more and more used in the context of nonlinearanalysis estimation Jauberthie et al. (2011) Raissi et al. 1. INTRODUCTION more and more used in the context of (2004) Raissi et are al. (2012) Drevelle and Bonnifait (2012) Reachability analysis is a classical problem in control interval nonlinear estimation Jauberthie et al. (2011) Raissi et al. Raissi et al. (2012) Drevelle and Bonnifait (2012) Reachability analysis is a classical problem in control (2004) nonlinear estimation Jauberthie et al. (2011) Raissi et al. dit Sandretto et al. (2014). theory Aubin and Frankowska (1990) Blanchini and MiRaissi et al. (2014). (2012) Drevelle and Bonnifait (2012) Reachability is a classical control dit Sandretto theory Aubin analysis and Frankowska (1990)problem Blanchiniin and Mi- (2004) (2004) Raissi et al. (2014). (2012) Drevelle and Bonnifait (2012) Reachability analysis a classical problem in and control ani (2007)Ramdani andis Nedialkov (2011)Bouissou etMial. dit Sandretto et al. theory Aubin and Frankowska (1990) Blanchini Assume now that there exists an unknown input vector ani (2007)Ramdani and Nedialkov (2011)Bouissou etMial. dit Sandretto et al. (2014). theory Aubin and Frankowska (1990) Blanchini (2014) and has several applications, for instanceand (i) to Assume now that there exists an unknown input vector ani (2007)Ramdani and Nedialkov (2011)Bouissou et al. u(k) ∈ [u] in the system which to a control (2014) and has several applications, for instance (i) to Assume now that therewhich existsmay an correspond unknown input vector ani (2007)Ramdani and applications, Nedialkov (2011)Bouissou et al. validate somehas properties of cyber-physic systems Konecny ∈ [u] in the system may to a control (2014) and several forsystems instance (i) to u(k) Assume now that there exists an correspond unknown input vector or a perturbation. We now have validate some properties of cyber-physic Konecny u(k) ∈ [u] in the system which may correspond to a control (2014) and has several applications, for instance (i) to et al. (2013)Taha and Duracz (2015), (ii) to ensure the or a perturbation. We now have validate some properties of cyber-physic systems Konecny u(k) ∈ [u] in the x(k system which may correspond to a control et al. (2013)Taha and Duracz (2015), (ii) to ensure the or a perturbation. We now have + 1) = f (x(k), u(k)). (5) validate some properties ofthe cyber-physic systems Konecny safeal.configuration during landing Bayen etensure al. (2007) et (2013)Taha and Duracz (2015),Bayen (ii) toet the or a perturbation. x(kWe + now 1) = have f (x(k), u(k)). (5) safe configuration during the landing al. (2007) et al.configuration (2013)Taha and Duracz (2015), (ii) toetensure the Thus the initial x(k + 1)set = fX(x(k), u(k)).also uncertain:(5) it feasible or (iii) to avoid collisions Desilles et al.Bayen (2012) with other 0 becomes safe during the landing al. (2007) becomes also uncertain: it Thus the initial feasible set X or (iii) to avoid collisions Desilles et al. (2012) with other x(k + 1) = f (x(k), u(k)). 0 safe configuration during the landing Bayen et with al. (2007) depends on the sequence u(k). Define the function ϕkk (5) as aircrafts. Reachability analysis allows us to guarantee that Thus becomes also uncertain: it the initial feasible set X or (iii) to avoid collisions Desilles et al. (2012) other 0 depends on the sequence u(k). Define the function ϕ as aircrafts. Reachability analysis allows us to guarantee that k it also uncertain: theon initial feasible set X0 becomes or toReachability avoid Desilles al. (2012) with aother follows. the(iii) system with acollisions givenanalysis control lawetwill always reach tar- Thus depends the sequence u(k). Define the function ϕ as aircrafts. allows us to guarantee that the systemReachability with a givenanalysis control law willusalways reach a that tar- follows.1 on the sequence u(k). Define the function ϕk as aircrafts. to guarantee get Aubin (2001). In this paper,allows we dealalways with the problem = f (x(0), u(0)) ϕ1 (x(0), u(0)) follows. the system(2001). with a given reach a tar- depends get Aubin In thiscontrol paper,law we will dealalways with the problem = f (x(0), u(0)) ϕ1 (x(0), u(0)) k+1 k the system with a given control law will reach a of computing the set X0 of all initial states xthe x(0),tarof follows. 0 = problem get Aubin (2001). In this paper, we deal with (x(0), u(0 : k)) = f (ϕ (x(0), u(0 : u(0)) k − 1)) , u(k)) ϕ (x(0), u(0)) f (x(0), ϕ k+11 k of computing the set X0 of all initial states xthe x(0), ¯of 0 = problem (x(0), u(0 : k)) = f (ϕ (x(0), u(0 : u(0)) k − 1)) , u(k)) ϕ ¯ get Aubin (2001). In this paper, we deal with (x(0), u(0)) = f (x(0), ϕ k+1 k aof discrete timethe system such that at time k¯ the state x(k) computing set X of all initial states x = x(0), of ¯ 0 0 (x(0), = f (ϕ (x(0), u(0 : k − 1)) , u(k)) ϕk+1u(0 aof discrete timethe system such that at time k¯ the state x(k) where : ku(0 − 1): k)) denotes the k sequence {u(0), . . . , u(k − computing set X of all initial states x = x(0), of ¯ is inside a given set Y. 0 0 where u(0 : k − 1) denotes the sequence {u(0), . . ,. u(k)) , u(k − (x(0), u(0 : k)) = f (ϕ (x(0), u(0 : k − 1)) ϕ If we ais discrete that at time k the state x(k) 1)}. inside a time givensystem set Y. such define ¯ where u(0 : k − 1) denotes the sequence {u(0), . . . , u(k − ais discrete time system such that at time k¯ the state x(k) 1)}. If we define inside a given set Y. where : k − 1) ¯denotes the {u(0),¯ . . . , u(k − ¯ sequence ¯ Weinside assume that set the system is described by k 1)}. If u(0 we is a given We assume that the Y. system is described by {x0define |∀u(0 : k−1) ∈ [u]k¯ , ϕkk¯ (x(0), u(0 : k−1)) ∈ Y} X⊂ ⊂ 0 = ¯ ¯ 1)}. If we define = {x |∀u(0 : k−1) ∈ [u] , ϕ (x(0), u(0 : k−1)) ∈ Y} X ¯ ¯ We assume that thex(k system is fdescribed 0 ¯k + 1) = (x(k)), by (1) 0 ¯k ⊂ k ⊃ k ¯ ¯ ¯ ¯ = {x |∀u(0 : k−1) ∈ [u] , ϕ (x(0), u(0 : k−1)) ∈ Y}, Y} X ¯ ¯ x(k + 1) = f (x(k)), (1) ¯ {x |∃u(0 : k−1) ∈ [u] , ϕ (x(0), u(0 : k−1)) ∈ 0 ¯ We assume thatn the system is described by n k 0 k 0 ⊃ k ⊂ k ¯ ¯ ¯ ¯ {x00|∃u(0 |∀u(0 :: k−1) k−1) ∈ ∈ [u] [u]k¯ ,,ϕ ϕk¯ (x(0), (x(0),u(0 u(0 :: k−1)) k−1)) ∈ ∈ Y}, Y} X⊃ + 1) = f (x(k)), (1) the state vector and f : Rn → Rnn is the where x(k) ∈ Rn is x(k 0 = {x ¯ ¯ (6) X0⊃ = {x0 |∃u(0 : k−1) ∈ [u]k¯ , ϕk¯ (x(0), u(0 : k−1)) ∈ Y}, the state vector and f : Rn → Rn is the where x(k) ∈ Rn is x(k + system. 1) = f (x(k)), (1) (6) ¯ ¯ evolution function of the If we define the function thethe state vector f : Rn the → Rfunction is the we where x(k)function ∈ Rn is of = {x0 |∃u(0 : k−1) ∈ [u] , ϕ (x(0), u(0 : k−1)) ∈ Y}, X0 have evolution system. If and we define (6) n the state f : R the → Rfunction is the where x(k)function ∈ R is of have evolution system. we (6) = f ◦ vector f ◦ · ·If · ◦and f , define (2) we ϕkkthe X⊂ X0 ⊂ X⊃ (7) ⊂ ⊃ we have 0 ⊂ 0 .. = system. f ◦ f ◦ · ·If · ◦we f , define the function (2) ϕkthe evolution function of X ⊂ X ⊂ X (7) 0 0 0 ⊂ ⊃ we have · · · ◦ f , (2) ϕk = f ◦ kf ◦times X ⊂ X ⊂ X . (7) 0 0 classically 0 · · · ◦ f , sets X⊂ and X⊃ are called the minimal and (2) The ϕ = f ◦ kf ◦times X⊂ ⊂ X0 ⊂ X⊃ . (7) 0 classically 0called The sets X0⊂ X0⊃ the minimal and we have k times 0 and 0 are ⊂ ⊃ ¯ the maximal backward reach set Jaulin and Walter (1994). we have −1 k times The sets X and X are classically called the minimal and k 0 0 ⊂ backward set Jaulin and Walter (1994). (3) the X0 = (ϕk¯ )−1 (Y) ⊃reach we have Themaximal sets andand X⊃ Xare classically called the minimal and Both setsX0X⊂ are set difficult toand compute, mainly (3) the X0 = (ϕk¯ )−1 (Y) maximal Jaulin Walter (1994). 0backward 0reach we have Both sets X⊂ and0 X⊃ are set difficult toand compute, mainly ) (Y) (3) which corresponds toXa0 = set(ϕinversion problem Jaulin and 0 0 ⊂ ⊃ ¯ the maximal backward reach Jaulin Walter (1994). due to the fact that existing methods cannot characterize k −1 Both sets X and X are difficult to compute, mainly which corresponds toXa0 = set(ϕinversion problem Jaulin and ) set (Y)inversion problem(3) 0 0 ⊂ that existing ⊃ due tosets the X fact methodsto cannot characterize Waltercorresponds (1993). Note is Both which to athat set the inversion problem problem Jaulin and X difficult compute, mainly fromtoinside the penumbra and Jaulin (2017) 0 and 0 are Desrochers the fact that existing methods cannot characterize Walter (1993). Note that the set inversion is due from inside the penumbra Desrochers and Jaulin (2017) which corresponds to a set inversion problem Jaulin and ⊃ ⊂ ∂ considered as much easier than the direct image problem Walter (1993). Note that the set inversion problem is due to the fact that existing methods cannot characterize which corresponds to the set X = X \X . Moreover, ⊃ ⊂ from inside the penumbra Desrochers and Jaulin (2017) ∂ considered as much easier than the direct image problem 0 0 0 which corresponds to the set X = X \X . Moreover, Walter (1993). Note that the set inversion problem is ⊂ ⊃ Mullier et al. (2013). Equivalently, we could have written 0 0 0 ⊃ ⊂ ∂or outside considered as (2013). much easier than thewe direct image problem from to inside the penumbra Desrochers and Jaulin (2017) even characterize inside X⊂ X.⊃ , existing which corresponds to the set = X \X Mullier et al. Equivalently, could have written 0X0 0 Moreover, 0 0 ⊃ ⊂ even to characterize inside X or outside X , existing ∂ considered as (2013). much easier than direct image problem which corresponds to the set ⊂ −1 thewe ¯ 0X 0 Moreover, ⊃ Mullier et Equivalently, could k = X \X . = f −1 = Y written (4) even X0al. 0 ¯ have 0 0 to characterize inside X or outside X , existing −1 (X1 ), X1 = f −1 (X2 ), . . . , Xk ¯ interval methods will need to 0⊂bisect inside the 0 ⊃ box [u]k Mullier et = (2013). f −1 (X1 ),Equivalently, X1 = f −1 (X2we ), . could . . , Xk¯ have = Y written (4) even X0al. interval methods will need to bisect inside the box [u] ¯ to characterize inside X or outside X , existing k 0 0This is not (X2 ), . . .to , Xk¯solve = Y a chain (4) to amounts control the accuracy of the characterization. and the Xcomputation of1 = X0f −1 0 = f −1 (X1 ), X interval methods will need tocharacterization. bisect inside theThis boxis [u] ¯ k to solve a chain to control the accuracy of the not and the Xcomputation of1 = X0f amounts = f (X ), X (X ), . . . , X = Y (4) ¯ 0 1 2 k the problem interval methods will todimension. bisect inside theThis boxisthere [u] of settheinversion problems. Inamounts both cases, satisfactory due to itsneed large Note that to control the accuracy of the characterization. not to solve a chain and computation of X 0 of problems. In both cases, the problem satisfactory to its of large dimension. Note This that isthere to control thedue accuracy the characterization. not to solve a chain andset theinversion computation of X can be solved using a Set Inversion Jaulin and satisfactory exists numerical methods Mitchell (2003)Note to compute an 0Inamounts due to its large dimension. that there of set problems. both approach cases, the problem can be inversion solved using a Set Inversion approach Jaulin and exists numerical methods Mitchell (2003)Note to compute an satisfactory due to its large dimension. that there of set inversion problems. In both cases, the problem Walter (1993) based on interval analysis. Set inversion and exists outer approximation of theMitchell maximal backward reach set numerical methods (2003) to compute an can be solved using a Set Inversion approach Jaulin Walter (1993) based on interval analysis. Set inversion and outer approximation of the maximal backward reach set numerical methods (2003) to compute an can be (1993) solved based using on a Set Inversion approach Jaulin and and exists outer approximation of theMitchell maximal backward reach set Walter interval analysis. Set inversion outer approximation of the maximal backward reach set Walter (1993) based on interval analysis. Set inversion and Copyright © 2017 IFAC 9674 ∗ ∗ ∗ ∗
Copyright © 2017, 2017 IFAC 9674Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © IFAC (International Federation of Automatic Control) Copyright © 2017 IFAC 9674 Peer review under responsibility of International Federation of Automatic Copyright © 2017 IFAC 9674Control. 10.1016/j.ifacol.2017.08.1297
Proceedings of the 20th IFAC World Congress Benoit Desrochers et al. / IFAC PapersOnLine 50-1 (2017) 9266–9271 Toulouse, France, July 9-14, 2017
X⊃ 0 , or an inner approximation of the minimal backward reach set X⊂ 0 , but the main difficulty in inside the penumbra, i.e., to compute accurately, an inner approximation ⊂ of X⊃ 0 and an outer approximation X0 . The main contribution of this paper is to show that by rewriting the problem as a chain of set inversion problem, we will be able to characterize, with an arbitrary accuracy ¯ k the feasible set X0 , without any bisection of [u] . An ⊂ inner approximation of the penumbra X∂0 = X⊃ 0 \X0 , which corresponds to the uncertainty of the approximation, will also be computed. This was not possible before except for some specific cases such as when the system is linear Girard (2003). The paper is organized as follows. Section 2 formulates the set inversion problem in the case where it depends on a parameter vector p which is assumed to be inside a box [p]. Section 3 introduces the specific case where the function to be inverted is a link function. This will allow us to build, in Section 4, a large class of set inversion problems that can be solved efficiently, without any bisection inside the parameter space. Test-cases related to reachability analysis are presented in Section 5. Section 6 concludes the paper. 2. PROBLEM Consider the set inversion problem (8) X = f −1 (Y) n m where f : R → R , and Y is a subset of Rm . Solving a set inversion problem consists of characterizing the solution set X from inside and from outside. We consider the uncertain case where f depends on a parameter vector p ∈ Rq . The function should thus be written as f (x, p), but for simplicity of notation when we will compose functions to form a chain, we will often write fp (x). In such a case, the solution set X also depends on the value of p and is denoted by X(p). We have: (9) X⊂ ⊂ X(p) ⊂ X⊃ , where fp−1 (Y) = {x|∀p ∈ [p], fp (x) ∈ Y} X⊂ =
and will be called a thick set inversion problem. To characterize X with an arbitrary accuracy, i.e., to characterize the sets X⊂ and X⊃ , existing interval methods, will need to bisect with respect to the p-space, which makes them inefficient. Moreover, they are not able to prove that a box is inside the penumbra. It has been show recently Desrochers and Jaulin (2017) how it can proved that a box is inside the penumbra. It will be shown in Section 3 that, in the specific case where the function fp (x) is a link function (i.e., for all x, and for any box [p], the set f[p] (x) is a box), bisections on [p] can be avoided. We now illustrate how this can be done on a simple example. Example 1. Consider the thick set inversion problem −1 X = f[p] ([y]) where fp (x) =
p∈[p]
Denote by (P(Rn ), ⊂) the powerset of Rn equipped with the inclusion ⊂ as an order relation. The powerset P(Rn ) is a complete lattice with respect to ⊂. The meet operator corresponds to the intersection and the join to the union. The pair [X⊂ , X⊃ ] is an interval in P(Rn ) which is called a thick set and denoted by X . A thick set partitions Rn into three zones: the clear zone X⊂ , the penumbra X∂ = X⊃ \X⊂ and the dark zone Rn \X⊃ . A thick set X is a sub-lattice of (P(Rn ), ⊂), i.e., if A ∈ X ,B ∈ X , then A ∩ B ∈ X and A ∪ B ∈ X .
Notation. The set inversion problem we want to solve will be written as −1 X = f[p] (Y) (10)
p1 p3
p2 p4
x1 x2
.
(11)
− − + + p1− x1 + p2− x2 , p1+ x1 + p2+ x2 −p3 x1 +−p4 x 2 , p 3 x+1 + p4 x+2 p1 x1 + p2 x2 p 1 x 1 + p2 x 2 = , − + p− x + p x p+ 1 2 1 + p4 x 2 3 3 x 4 − + = f[p] (x), f[p] (x)
f (x, [p]) =
which corresponds to a box. If x ∈ [x] = [0, 1] × [2, 3], f − (x, [p]) belongs to the box 2 · [0, 1] + 2 · [2, 3] − ([x]) = f[p] 4 · [0, 1] − 6 · [2, 3] [0, 2] + [4, 6] = [0, 4] + [−18, −12] [4, 8] = [−18, −8]
fp−1 (Y) = {x|∃p ∈ [p], fp (x) ∈ Y} .
The set X∂ = X⊃ \X⊂ is called the penumbra Grimson and Lozano-P´erez (1985) and contains the unknown boundary of X.
Assume that [p] = [2, 3] × [2, 5] × [4, 5] × [−6, −1] and [y] = [5, 19] × [−7, 11]. Take a box [x] = [0, 1] × [2, 3] and let us show how is can be proved that it is inside the penumbra. Note that for simplicity, we have taken all intervals with a constant sign, but the method can be made much more general if we use modal interval analysis Goldsztejn et al. (2005) Vinas et al. (2006) or symbolic interval arithmetic Jaulin and Chabert (2010). Inside [x], we have
p∈[p]
X⊃ =
9267
(12)
and f + (x, [p]) belongs to the box 3 · [0, 1] + 5 · [2, 3] + ([x]) = f[p] 5 · [0, 1] − 1 · [2, 3] [0, 3] + [10, 15] = [0, 5] + [−3, −2] [10, 18] = . [−3, 3]
(13)
− + ([x]) ∩ [y] = ∅ and f[p] ([x]) ⊂ [y], we Since f[p] conclude that [x] is inside the penumbra. This is illustrated by Figure 1.
9675
Proceedings of the 20th IFAC World Congress 9268 Benoit Desrochers et al. / IFAC PapersOnLine 50-1 (2017) 9266–9271 Toulouse, France, July 9-14, 2017
0 -0.5
0
0.5
1
1.5
8
3.4
6
3.2
5
10
15
20
25
10
3.6
4 2
3 0 2.8 2.6 2.4 2.2 2
-2 -4 -6 -8 -10 -12
1.8 1.6 1.4
-14 -16 -18
+ Fig. 1. For [x] = [0, 1]×[2, 3] the set f[p] (blue dots) is inside − [y] whereas f[p] (red dots) is outside. We conclude that [x] is inside the penumbra
Fig. 2. Left: Classical interval methods accumulate on the yellow thick boundary (the penumbra). Right: the method we propose here allows a fast treatment of the penumbra. The frame box is [−5, 11]2 3. LINK FUNCTIONS
In the specific case where the function fp (x) is a link function (see Section 3 for more details), for all x, the set − + f[p] (x) is a box, denoted by f[p] (x), f[p] (x) , with lower
− + bound f[p] (x) ∈ Rm and upper bound f[p] (x) ∈ Rm . These two bounds can be described by an algorithm and two inclusion functions can be obtained for them using the rules of interval computation Kearfott (1996). This means that if x ∈ [x], we can obtain, under the form of an − + algorithm, two boxes f[p] ([x]) and f[p] ([x]) such that
− − + + f[p] ([x]) and f[p] ([x]). ([x]) ∈ f[p] ([x]) ∈ f[p] To characterize the thick set X of the thick set inversion −1 problem X = f[p] (Y), we start with a huge box (large enough to contain X⊃ ) and we cut it into several other boxes all stored inside a list L. For each box [x] of L, and in − + we compute the two boxes f[p] ([x]) and f[p] ([x]). From
these two boxes, we are able to test if [x] ⊂ X⊂ or if [x] ∩ X⊃ = ∅ or if [x] is included inside the penumbra. If nothing can be concluded, the box [x] is bisected and the two resulting boxes are stored inside the list. The treatment is applied to all boxes of the list until all boxes become too small (i.e., with a width smaller than a threshold ε). Note that, if we are able to conclude that a box [x] is inside the penumbra, no bisection would occur.
For Example 1, the method generates Figure 2 (right) which is an approximation of a thick set with the inner part (red), the outer part (blue) and the penumbra (orange). A classical interval method would not able to characterize the penumbra and yields Figure 2 (left). Since the resulting algorithm bisects in the x-sp ace, its complexity is exponential with respect to the dimension of x. Its application is thus limited to low dimensions. If do not bisect on the p-space (as we suggest in this paper using chains and link functions), the complexity with respect to p can be considered as polynomial.
This section considers the case of thick set-inversion problems in the case where the functions to be inverted have some good properties. These functions will be the links to build more complex functions to be inverted. Definition 2. Link function. If, for a box [p] ⊂ Rq , the set f[p] (x) = f (x, [p]) = {y ∈ Rp | ∃p ∈ [p] , y = fp (x)} is a box, then the function f is said to be a link. Links will be composed later in Section 4 in order to build a chain. If m = 1 then the function f (x, p) becomes scalar. Note that all f (x, p) that are scalar and continuous with respect to p, are links. Due to their specific box-shaped structure, link functions can be inverted with respect to x without bisections with respect to the uncertain parameter box [p]. Example 3. The function f (x, p) = 20e−x1 p − 8e−x2 p (14) is a scalar link function, due to the continuity of f with respect to p. Now, to get a computable expression for f (x, [p]), we need to compute the extremum of the function inside the interval [p] = [p− , p+ ]. Since ∂f (x, p¯) =0 ∂p ⇔ −20x1 · e−x1 p¯ + 8x2 · e−x2 p¯ = 0 2x2 (15) ⇔ e(x2 −x1 )p¯ = 5x 1 1 2x2 , ⇔ p¯ = ln x2 − x1 5x1 inside the interval [p] = [p− , p+ ], f (x, p), may have the extremum y¯ = f (x, p¯). As a result, if we define y − = f (x, p− ), y + = f (x, p+ ), we have − + min{y , y , y¯}, max{y − , y + , y¯} if p¯ ∈ [p] f (x, [p]) = if p¯ ∈ / [p] min{y − , y + }, max{y − , y + } Equivalently, − + (x), f[p] (x) (16) f (x, [p]) = f[p] where
9676
− (x) f[p]
=
min{y − , y + } if p¯ ∈ / [p] min{y − , y + , y¯} if p¯ ∈ [p]
(17)
Proceedings of the 20th IFAC World Congress Benoit Desrochers et al. / IFAC PapersOnLine 50-1 (2017) 9266–9271 Toulouse, France, July 9-14, 2017
and + (x) f[p]
=
/ [p] max{y − , y + } if p¯ ∈ max{y − , y + , y¯} if p¯ ∈ [p]
9269
(18)
− + (x) and f[p] (x) belong to R. If x ∈ [x], Note that both f[p] we have − + p] ∩ [p] = ∅ min{ − − −y , +y } if [¯ ([x]) = (x) ∈ f[p] f[p] y ]} otherwise min{ y , y , [¯
and
− + p] ∩ [p] = ∅ max{ −y , +y } if [¯ y ]} otherwise Fig. 3. Left: measurement boxes [t ] × [y ]. Right: Approxmax{ y , y , [¯ i i where imation of the thick set Q in the parameter space 1 2 [x2 ] ln [¯ p] = Example 7. Consider the problem of estimating the pa[x2 ] − [x1 ] 5 [x1 ] y] = f ([x] , [¯ p] ∩ [p]) (19) rameters q1 and q2 of the model [¯ − − y(q, t) = 20e−q1 t − 8e−q2 t . f ([x] , p ) y + = y = f ([x] , p+ ). We assume that 10 measurements yi have been collected at Definition 4. Vector components. The family of vectors time ti . The uncertainties on the pair (ti , yi ) is represented {p1 , p2 , . . . , p } correspond to vector components of the by intervals as represented by the table below and Figure vector p if each pi is a subvector of p and the indexes are 3 (left). We are interested by the thick set all disjoint. For instance if p = (p1 , . . . , p9 ), then {p1 , p2 }, Q = q ∈ R2 |∀i, y(q, ti ) ∈ [yi ] . where p1 = (p2 , p9 ) and p2 = (p8 , p1 , p3 ), correspond to If we set vector components of p. Note that the two index sets {2, 9} y(q, [t1 ]) and {8, 1, 3} are disjoint. .. f[t] (q) = . The fact that index sets are disjoint implies in independeny(q, [t10 ]) cies between parameters. The conservatism related to the [y] = [y10 ] × · · · × [y10 ] dependency effect, well known in interval analysis Moore [t] = [t1 ] × · · · × [t10 ] (1979), is thus avoided. In practice, this will allow us to then, we have avoid bisection in the {p1 , p2 , . . . , p }-space. −1 Q = f[t] ([y]), Proposition 5. Consider scalar link functions, fi (x, pi ) , where vectors {p1 , p2 , . . . , p } are vector components of which is thick set inversion problem with the link function the vector p. The function f[t] (q) to be inverted. The inversion yields the approxima tion of Figure 3 (right). f1 (x, p1 ) . .. (20) f (x, p) = i [ti ] [yi ] 1 [0.25, 1.25] [2.7,12.1] f (x, p ) 2 [1, 2] [1.04,7.14] is a link. 3 [1.75, 2.75] [-0.13,3.61] 4 [2.5, 3.5] [-0.95,1.15] Proof. Consider the box [p] and its box components 5 [5.5, 6.5] [-4.85,-0.29] [p1 ] , . . . , [p ]. The set f (x, [p]) is the Cartesian product 6 [8.5, 9.5] [-5.06,-0.36] of intervals: 7 [12.5, 13.5] [-4.1,-0.04] f (x, [p]) = f1 (x, [p1 ]) × · · · × f (x, [p ]) (21) 8 [16.5, 17.5] [-3.16,0.3] 9 [20.5, 21.5] [-2.5,0.51] which corresponds to a box. 10 [24.5, 25.5] [-2,0.6] Example 6. The function p1 x1 x2 + p1 p2 sin(x2 ) 4. CHAIN f (x, p) = (22) p3 x1 + p3 p4 cos(x1 )x2 + f[p] (x)
+ ∈ f[p] ([x]) =
is a link. Now, if we replace p3 by p2 , a dependency occurs between the two components of f and the function is not a link anymore.
Consequence. If fp (x) is a link, for all x, the set f[p] (x) is a box. If x ∈ [x], the lower and upper bounds of the − box f[p] (x) is included inside the boxes f[p] ([x]) and + f[p] ([x]). As a consequence, the methodology explained in Section 2 can be applied to find an inner and outer −1 approximations of the thick set X = f[p] (Y), with an approximation of the penumbra.
Definition 8. Chain. The function ϕ(x, p) = ϕp (x) is said to be a chain if it can be written as (23) ϕp (x) = fp ◦ · · · ◦ fp22 ◦ fp11 (x)
where {p1 , p2 , . . . , p } are vector components of p and the fpkk are links.
Note that the parameter vectors p1 , p2 , . . . , p are independent. This implies that the pessimism due to the dependency effect will not take place during the resolution. Moreover, the fact that fpkk are links implies that pessimism due to the wrapping effect (of a set by a box) will not exist.
9677
Proceedings of the 20th IFAC World Congress 9270 Benoit Desrochers et al. / IFAC PapersOnLine 50-1 (2017) 9266–9271 Toulouse, France, July 9-14, 2017
Example 9. The function
sin(p3 x21 ) + p4 x ϕp (x) = p5 (p1 e 2 + p1 p2 x1 x2 ) p3 x21
(24)
is not a link. Indeed, p3 occurs on both components of ϕp . Now, if we define p 1 e x 2 + p 1 p2 x 1 x 2 1 fp1 (x) = (25) p3 x21 fp22 (z)
=
sin(z2 ) + p4 p 5 z1 z2
(26)
with p1 = (p1 , p2 , p3 ) and p2 = (p4 , p5 ) then, we have ϕp (x) = fp22 ◦ fp11 (x).
(27)
Since both fp22 , fp11 are links, ϕp (x) is a chain. Consequence. Solving a thick set inversion problem X = ϕ−1 [p] (Y) in the case where ϕ is a chain can be done by solving several thick set inversion problems, without any bisection inside the box [p]. Rewriting a chain as a composition of link functions can probably be done by symbolic methods, but we do not know any approach or algorithm to perform this task. Of course, we should not allow the dimension of the intermediate spaces to be arbitrary large compare to the dimension of x if we want to be efficient. The following section provides some test-cases of the thick inversion of chains, applied to reachability problems.
5. TEST-CASES
Fig. 4. Approximation of the sets X(0), X(1), X(2), X(3) Test-case 2. Consider the nonlinear system x (k) + x2 (k) · u (k) 1 1 2 x1 (k + 1) = 1 (31) x2 (k + 1) · x1 (k) · x2 (k) + u2 (k) 2 with u1 ∈ [1, 2] and u2 ∈ [−2, −1]. The set Y is assumed to be a centred disk which has to be reached at time k¯ = 3. Solving the chain of thick set inversion problem generates Figure 5. The frame boxes are [−12, 10] × [−12, 10] for all X (i).
Test-case 1. Consider the linear system x(k + 1) = A · x(k) where A(k) ∈
[2.5, 3] [2, 3] [4, 4.5] [−3, −2]
(28)
(29)
plays the role of the parameter vector p. Assume that we want to reach the target Y = [4, 20] × [−8, 12] (30) ¯ at time k = 3. Solving the chain of thick set inversion problems yields the approximation of the thick sets X (0), X (1), X (2), X (3) as depicted on Figure 4. Note that X (3) corresponds to the box Y. We first compute the thick set X (2) = f −1 [A] ( X (3)). Note that X (3) is here a thin set (i.e., a classical set) whereas X (2) is thick (i.e., with a penumbra). Then we compute X (1) = f −1 [A] ( X (2)) and finally, we compute X (0) = −1 f [A] ( X (1)). At each step, the penumbra inflates, due to the accumulation of uncertainties. The frame boxes are [−0.6, 0.8] × [−0.4, 1] for X (0), [−1.5, 3]×[−2.5, 2.2] for X (1),[−4, 8]×[−3, 10] for X (2), and [−5, 27] × [−15, 20] for X (3).
Fig. 5. Approximation of the sets X(0), X(1), X(2), X(3) 6. CONCLUSION In this paper, we have presented a new interval based approach to compute the set X0 of initial states that will
9678
Proceedings of the 20th IFAC World Congress Benoit Desrochers et al. / IFAC PapersOnLine 50-1 (2017) 9266–9271 Toulouse, France, July 9-14, 2017
reach a target Y in a finite time. The dynamic system that is considered is discrete time and uncertain (the uncertainties are represented by intervals). The principle is to transform the problem into a chain of set inversion problems composed with links, i.e., a function with nice properties for set inversion. Among the problems that remain to be solved, we need (1) to find symbolic methods for the decomposition of a chain into links with low dimensions; (2) to extend the approach to continuous systems described by differential equations and to study how the penumbra can be characterized from inside in such a case; and (3) to solve the chain of set inversion problems in parallel and not sequentially as done here. REFERENCES Aubin, J.P. (2001). Viability Kernels and Capture Basins of Sets Under Differential Inclusions. SIAM Journal on Control and Optimization, 40(3), 853–881. doi: 10.1137/S036301290036968X. Aubin, J. and Frankowska, H. (1990). Set-Valued Analysis. Birkh¨ auser, Boston. Bayen, A.M., Mitchell, I.M., Osihi, M.K., and Tomlin, C.J. (2007). Aircraft Autolander Safety Analysis Through Optimal Control-Based Reach Set Computation. Journal of Guidance, Control, and Dynamics, 30(1), 68–77. doi:10.2514/1.21562. Blanchini, F. and Miani, S. (2007). Set-Theoretic Methods in Control. Springer Science & Business Media. Bouissou, O., Chapoutot, A., Djaballah, A., and Kieffer, M. (2014). Computation of parametric barrier functions for dynamical systems using interval analysis. In IEEE CDC. Los Angeles, United States. Desilles, A., Zidani, H., and Cruck, E. (2012). Collision analysis for an UAV. In AIAA Guidance, Navigation, and Control Conference, 13–16. American Institute of Aeronautics and Astronautics. Desrochers, B. and Jaulin, L. (2017). Computing a guaranteed approximation the zone explored by a robot. IEEE Transaction on Automatic Control, 62(1), 425– 430. dit Sandretto, J.A., Trombettoni, G., Daney, D., and Chabert, G. (2014). Certified calibration of a cabledriven robot using interval contractor programming. In F. Thomas and A.P. Gracia (eds.), Computational Kinematics, Mechanisms and Machine Science, Spinger. Drevelle, V. and Bonnifait, P. (2012). igps: Global positioning in urban canyons with road surface maps. IEEE Intelligent Transportation Systems Magazine, 4(3), 6– 18. Girard, A. (2003). Computation and stability analysis of limit cycles in piecewise linear hybrid systems. In 1st IFAC Conference on Analysis and Design of Hybrid Systems, 181–186. Goldsztejn, A., Daney, D., Rueher, M., and Taillibert, P. (2005). Modal intervals revisited: a mean-value extension to generalized intervals. In Proceedings of the Quantification in Constraint Programming, Barcelona, Spain, 2005 (QCP 2005). Grimson, W.E.L. and Lozano-P´erez, T. (1985). Recognition and localization of overlapping parts in two and three dimensions. In IEEE International Conference on Robotics and Automation, 61–66.
9271
Jauberthie, C., Verdi`ere, N., and Trav´e-Massuy`es, L. (2011). Set-membership identifiability: definitions and analysis. In In Proceedings of the 18th IFAC World Congress, Milan, Italie, volume 00, 12024–12029. aaa. Jaulin, L. and Chabert, G. (2010). Resolution of nonlinear interval problems using symbolic interval arithmetic. Engineering Applications of Artificial Intelligence, 23(6), 1035–1049. Jaulin, L. and Walter, E. (1993). Guaranteed nonlinear parameter estimation via interval computations. Interval Computation, 61–75. Jaulin, L. and Walter, E. (1994). Modeling Techniques for Uncertain Systems, chapter Set Inversion, with Application to Guaranteed Nonlinear Estimation and Robust Control, 3–20. Birkh¨auser. Kearfott, R.B. (1996). Rigorous Global Search: Continuous Problems. Kluwer, Dordrecht, the Netherlands. Konecny, M., Taha, W., Duracz, J., Duracz, A., and Ames, A. (2013). Enclosing the behavior of a hybrid system up to and beyond a zeno point. In Cyber-Physical Systems, Networks, and Applications (CPSNA). Mitchell, I.M. (2003). Application of level set methods to control and reachability problems in continuous and hybrid systems. Moore, R.E. (1979). Methods and Applications of Interval Analysis. SIAM, Philadelphia, PA. Mullier, O., Goubault, E., Kieffer, M., and Putot, S. (2013). General inner approximation of vector-valued functions. Reliable Computing, 18(0), 117–143. Raissi, T., Efimov, D., and Zolghadri, A. (2012). Interval state estimation for a class of nonlinear systems. IEEE Transactions on Automatic Control, 57(1), 260–265. Raissi, T., Ramdani, N., and Candau, Y. (2004). Set membership state and parameter estimation for systems described by nonlinear differential equations. Automatica, 40, 1771–1777. Ramdani, N. and Nedialkov, N. (2011). Computing Reachable Sets for Uncertain Nonlinear Hybrid Systems using Interval Constraint Propagation Techniques. Nonlinear Analysis: Hybrid Systems, 5(2), 149–162. Taha, W. and Duracz, A. (2015). Acumen: An opensource testbed for cyber-physical systems research. In CYCLONE’15. Vinas, P.H., Sainz, M.A., Vehi, J., and Jaulin, L. (2006). Quantified set inversion algorithm with applications to control. Reliable computing, 11(5), 369–382.
9679
ScienceDirect
IFAC PapersOnLine 50-1 Inversion (2017) 9266–9271 Problems Chain of Set Chain of Set Inversion Problems Chain of Set Inversion Problems Application to reachability analysis Chain of Set Inversion Problems Application to reachability Application to reachability analysis analysis Application to reachability analysis ∗ ∗∗ Benoit Desrochers Luc Jaulin ∗ ∗∗
Benoit Desrochers ∗ Luc Jaulin ∗∗ Benoit Desrochers ∗ Luc Jaulin ∗∗ Benoit DesrochersTechniques Luc Jaulin Direction G´en´erale de l’Armement, Navales Brest, BCRM Direction G´en´erale de l’Armement, Techniques Navales Brest, BCRM Brest, 29240 BrestTechniques Cedex 9, (e-mail: Direction G´en´ede rale de l’Armement, Navales Brest, BCRM Brest, 29240 BrestTechniques Cedex 9, (e-mail: Direction G´en´ede rale de l’Armement, Navales Brest, BCRM [email protected]) de Brest, 29240 Brest Cedex 9, (e-mail: [email protected]) ∗∗ de Brest, 29240 Brest Cedex 9, (e-mail: ENSTA-Bretagne, Lab-STICC, 2 rue Fran¸ c ois Verny, 29806 Brest, [email protected]) ∗∗ Lab-STICC, 2 rue Fran¸cois Verny, 29806 Brest, ∗∗ ENSTA-Bretagne, [email protected]) (e-mail: [email protected] )Verny, 29806 Brest, ENSTA-Bretagne, Lab-STICC, 2 rue Fran¸ c ois (e-mail: [email protected] ∗∗ ENSTA-Bretagne, Lab-STICC, 2 rue Fran¸cois))Verny, 29806 Brest, (e-mail: [email protected] (e-mail: [email protected] ) Abstract: This paper deals with the set inversion problem X = f −1 (Y) in the case where Abstract: This paper deals with the set inversion qproblem X = f −1 (Y) in the case where m −1 to f : Rnn → RThis on a parameter vector p ∈ Rqproblem which isX known a box [p]. (Y)beininside the case where Abstract: paper deals with the set inversion = f m depends f : Rn → RThis depends on a parameter vector p ∈ Rqproblem which isX known −1 to be inside a box [p]. m (Y) ininside the case where Abstract: paper deals with the set inversion = f We show that for a large class of problems, wep can obtain anis accurate approximation of the fWe : R → R depends on a parameter vector ∈ R which known to be a box show that for a large of problems, can obtain anis accurate of [p]. the n q f : R → Rm depends on class a parameter vectorwe ∈R which known toapproximation be a box [p]. solution set, without bisecting inofthe p-space. Top do this, symbolic methods areinside required to cast We show that for a large class problems, we can obtain an accurate approximation of the solution set, without bisecting in the p-space. To do this, symbolic methods are required to cast We show that for a into large can obtain an accurate approximation of cast the our initial problem a class chain ofproblems, set-inversion problems, the links of which have some nice solution set, without bisecting inofthe p-space. we To do this, symbolic methods are required to our initial problem into a chain of set-inversion problems, the links of which have some nice solution set, without bisecting in an the p-space. To we do this, symbolic methods required to cast properties with respect toa p. As application, consider thelinks problem ofare computing the set our initial problem into chain of set-inversion problems, the of which have some nice properties with respect toa p. As an application, we considerthe thelinks problem of computing thenice set our problem chain of set-inversion problems, have of allinitial initialwith states ofinto an uncertain discrete-time state system that reachofa which target set Y some in athe given properties respect to p. As an application, we consider the problem of computing set of all initial states of an uncertain discrete-time state system that reach a target set Y in a given properties with respect to p. As an application,state we consider the problem of computing the set time. of all initial states of an uncertain discrete-time system that reach a target set Y in a given time. of all initial states of an uncertain discrete-time state system that reach a target set Y in a given time. © 2017, IFAC (Internationalanalysis, Federation of Automatic Control) Hostingsystems; by Elsevier Ltd. All error rights reserved. time. Keywords: Reachability Interval analysis, Quantized Bounded Keywords: Reachability analysis, Interval analysis, Quantized systems; Bounded error identification, Keywords: Reachability analysis, Interval analysis, Quantized systems; Bounded error identification, Keywords: Reachability analysis, Interval analysis, Quantized systems; Bounded error identification, identification, 1. INTRODUCTION interval analysis are more and more used in the context of 1. INTRODUCTION interval analysis are more and more used in the context of nonlinearanalysis estimation Jauberthie et al. (2011) Raissi et al. 1. INTRODUCTION interval are more and more used in the context of nonlinearanalysis estimation Jauberthie et al. (2011) Raissi et al. 1. INTRODUCTION more and more used in the context of (2004) Raissi et are al. (2012) Drevelle and Bonnifait (2012) Reachability analysis is a classical problem in control interval nonlinear estimation Jauberthie et al. (2011) Raissi et al. Raissi et al. (2012) Drevelle and Bonnifait (2012) Reachability analysis is a classical problem in control (2004) nonlinear estimation Jauberthie et al. (2011) Raissi et al. dit Sandretto et al. (2014). theory Aubin and Frankowska (1990) Blanchini and MiRaissi et al. (2014). (2012) Drevelle and Bonnifait (2012) Reachability is a classical control dit Sandretto theory Aubin analysis and Frankowska (1990)problem Blanchiniin and Mi- (2004) (2004) Raissi et al. (2014). (2012) Drevelle and Bonnifait (2012) Reachability analysis a classical problem in and control ani (2007)Ramdani andis Nedialkov (2011)Bouissou etMial. dit Sandretto et al. theory Aubin and Frankowska (1990) Blanchini Assume now that there exists an unknown input vector ani (2007)Ramdani and Nedialkov (2011)Bouissou etMial. dit Sandretto et al. (2014). theory Aubin and Frankowska (1990) Blanchini (2014) and has several applications, for instanceand (i) to Assume now that there exists an unknown input vector ani (2007)Ramdani and Nedialkov (2011)Bouissou et al. u(k) ∈ [u] in the system which to a control (2014) and has several applications, for instance (i) to Assume now that therewhich existsmay an correspond unknown input vector ani (2007)Ramdani and applications, Nedialkov (2011)Bouissou et al. validate somehas properties of cyber-physic systems Konecny ∈ [u] in the system may to a control (2014) and several forsystems instance (i) to u(k) Assume now that there exists an correspond unknown input vector or a perturbation. We now have validate some properties of cyber-physic Konecny u(k) ∈ [u] in the system which may correspond to a control (2014) and has several applications, for instance (i) to et al. (2013)Taha and Duracz (2015), (ii) to ensure the or a perturbation. We now have validate some properties of cyber-physic systems Konecny u(k) ∈ [u] in the x(k system which may correspond to a control et al. (2013)Taha and Duracz (2015), (ii) to ensure the or a perturbation. We now have + 1) = f (x(k), u(k)). (5) validate some properties ofthe cyber-physic systems Konecny safeal.configuration during landing Bayen etensure al. (2007) et (2013)Taha and Duracz (2015),Bayen (ii) toet the or a perturbation. x(kWe + now 1) = have f (x(k), u(k)). (5) safe configuration during the landing al. (2007) et al.configuration (2013)Taha and Duracz (2015), (ii) toetensure the Thus the initial x(k + 1)set = fX(x(k), u(k)).also uncertain:(5) it feasible or (iii) to avoid collisions Desilles et al.Bayen (2012) with other 0 becomes safe during the landing al. (2007) becomes also uncertain: it Thus the initial feasible set X or (iii) to avoid collisions Desilles et al. (2012) with other x(k + 1) = f (x(k), u(k)). 0 safe configuration during the landing Bayen et with al. (2007) depends on the sequence u(k). Define the function ϕkk (5) as aircrafts. Reachability analysis allows us to guarantee that Thus becomes also uncertain: it the initial feasible set X or (iii) to avoid collisions Desilles et al. (2012) other 0 depends on the sequence u(k). Define the function ϕ as aircrafts. Reachability analysis allows us to guarantee that k it also uncertain: theon initial feasible set X0 becomes or toReachability avoid Desilles al. (2012) with aother follows. the(iii) system with acollisions givenanalysis control lawetwill always reach tar- Thus depends the sequence u(k). Define the function ϕ as aircrafts. allows us to guarantee that the systemReachability with a givenanalysis control law willusalways reach a that tar- follows.1 on the sequence u(k). Define the function ϕk as aircrafts. to guarantee get Aubin (2001). In this paper,allows we dealalways with the problem = f (x(0), u(0)) ϕ1 (x(0), u(0)) follows. the system(2001). with a given reach a tar- depends get Aubin In thiscontrol paper,law we will dealalways with the problem = f (x(0), u(0)) ϕ1 (x(0), u(0)) k+1 k the system with a given control law will reach a of computing the set X0 of all initial states xthe x(0),tarof follows. 0 = problem get Aubin (2001). In this paper, we deal with (x(0), u(0 : k)) = f (ϕ (x(0), u(0 : u(0)) k − 1)) , u(k)) ϕ (x(0), u(0)) f (x(0), ϕ k+11 k of computing the set X0 of all initial states xthe x(0), ¯of 0 = problem (x(0), u(0 : k)) = f (ϕ (x(0), u(0 : u(0)) k − 1)) , u(k)) ϕ ¯ get Aubin (2001). In this paper, we deal with (x(0), u(0)) = f (x(0), ϕ k+1 k aof discrete timethe system such that at time k¯ the state x(k) computing set X of all initial states x = x(0), of ¯ 0 0 (x(0), = f (ϕ (x(0), u(0 : k − 1)) , u(k)) ϕk+1u(0 aof discrete timethe system such that at time k¯ the state x(k) where : ku(0 − 1): k)) denotes the k sequence {u(0), . . . , u(k − computing set X of all initial states x = x(0), of ¯ is inside a given set Y. 0 0 where u(0 : k − 1) denotes the sequence {u(0), . . ,. u(k)) , u(k − (x(0), u(0 : k)) = f (ϕ (x(0), u(0 : k − 1)) ϕ If we ais discrete that at time k the state x(k) 1)}. inside a time givensystem set Y. such define ¯ where u(0 : k − 1) denotes the sequence {u(0), . . . , u(k − ais discrete time system such that at time k¯ the state x(k) 1)}. If we define inside a given set Y. where : k − 1) ¯denotes the {u(0),¯ . . . , u(k − ¯ sequence ¯ Weinside assume that set the system is described by k 1)}. If u(0 we is a given We assume that the Y. system is described by {x0define |∀u(0 : k−1) ∈ [u]k¯ , ϕkk¯ (x(0), u(0 : k−1)) ∈ Y} X⊂ ⊂ 0 = ¯ ¯ 1)}. If we define = {x |∀u(0 : k−1) ∈ [u] , ϕ (x(0), u(0 : k−1)) ∈ Y} X ¯ ¯ We assume that thex(k system is fdescribed 0 ¯k + 1) = (x(k)), by (1) 0 ¯k ⊂ k ⊃ k ¯ ¯ ¯ ¯ = {x |∀u(0 : k−1) ∈ [u] , ϕ (x(0), u(0 : k−1)) ∈ Y}, Y} X ¯ ¯ x(k + 1) = f (x(k)), (1) ¯ {x |∃u(0 : k−1) ∈ [u] , ϕ (x(0), u(0 : k−1)) ∈ 0 ¯ We assume thatn the system is described by n k 0 k 0 ⊃ k ⊂ k ¯ ¯ ¯ ¯ {x00|∃u(0 |∀u(0 :: k−1) k−1) ∈ ∈ [u] [u]k¯ ,,ϕ ϕk¯ (x(0), (x(0),u(0 u(0 :: k−1)) k−1)) ∈ ∈ Y}, Y} X⊃ + 1) = f (x(k)), (1) the state vector and f : Rn → Rnn is the where x(k) ∈ Rn is x(k 0 = {x ¯ ¯ (6) X0⊃ = {x0 |∃u(0 : k−1) ∈ [u]k¯ , ϕk¯ (x(0), u(0 : k−1)) ∈ Y}, the state vector and f : Rn → Rn is the where x(k) ∈ Rn is x(k + system. 1) = f (x(k)), (1) (6) ¯ ¯ evolution function of the If we define the function thethe state vector f : Rn the → Rfunction is the we where x(k)function ∈ Rn is of = {x0 |∃u(0 : k−1) ∈ [u] , ϕ (x(0), u(0 : k−1)) ∈ Y}, X0 have evolution system. If and we define (6) n the state f : R the → Rfunction is the where x(k)function ∈ R is of have evolution system. we (6) = f ◦ vector f ◦ · ·If · ◦and f , define (2) we ϕkkthe X⊂ X0 ⊂ X⊃ (7) ⊂ ⊃ we have 0 ⊂ 0 .. = system. f ◦ f ◦ · ·If · ◦we f , define the function (2) ϕkthe evolution function of X ⊂ X ⊂ X (7) 0 0 0 ⊂ ⊃ we have · · · ◦ f , (2) ϕk = f ◦ kf ◦times X ⊂ X ⊂ X . (7) 0 0 classically 0 · · · ◦ f , sets X⊂ and X⊃ are called the minimal and (2) The ϕ = f ◦ kf ◦times X⊂ ⊂ X0 ⊂ X⊃ . (7) 0 classically 0called The sets X0⊂ X0⊃ the minimal and we have k times 0 and 0 are ⊂ ⊃ ¯ the maximal backward reach set Jaulin and Walter (1994). we have −1 k times The sets X and X are classically called the minimal and k 0 0 ⊂ backward set Jaulin and Walter (1994). (3) the X0 = (ϕk¯ )−1 (Y) ⊃reach we have Themaximal sets andand X⊃ Xare classically called the minimal and Both setsX0X⊂ are set difficult toand compute, mainly (3) the X0 = (ϕk¯ )−1 (Y) maximal Jaulin Walter (1994). 0backward 0reach we have Both sets X⊂ and0 X⊃ are set difficult toand compute, mainly ) (Y) (3) which corresponds toXa0 = set(ϕinversion problem Jaulin and 0 0 ⊂ ⊃ ¯ the maximal backward reach Jaulin Walter (1994). due to the fact that existing methods cannot characterize k −1 Both sets X and X are difficult to compute, mainly which corresponds toXa0 = set(ϕinversion problem Jaulin and ) set (Y)inversion problem(3) 0 0 ⊂ that existing ⊃ due tosets the X fact methodsto cannot characterize Waltercorresponds (1993). Note is Both which to athat set the inversion problem problem Jaulin and X difficult compute, mainly fromtoinside the penumbra and Jaulin (2017) 0 and 0 are Desrochers the fact that existing methods cannot characterize Walter (1993). Note that the set inversion is due from inside the penumbra Desrochers and Jaulin (2017) which corresponds to a set inversion problem Jaulin and ⊃ ⊂ ∂ considered as much easier than the direct image problem Walter (1993). Note that the set inversion problem is due to the fact that existing methods cannot characterize which corresponds to the set X = X \X . Moreover, ⊃ ⊂ from inside the penumbra Desrochers and Jaulin (2017) ∂ considered as much easier than the direct image problem 0 0 0 which corresponds to the set X = X \X . Moreover, Walter (1993). Note that the set inversion problem is ⊂ ⊃ Mullier et al. (2013). Equivalently, we could have written 0 0 0 ⊃ ⊂ ∂or outside considered as (2013). much easier than thewe direct image problem from to inside the penumbra Desrochers and Jaulin (2017) even characterize inside X⊂ X.⊃ , existing which corresponds to the set = X \X Mullier et al. Equivalently, could have written 0X0 0 Moreover, 0 0 ⊃ ⊂ even to characterize inside X or outside X , existing ∂ considered as (2013). much easier than direct image problem which corresponds to the set ⊂ −1 thewe ¯ 0X 0 Moreover, ⊃ Mullier et Equivalently, could k = X \X . = f −1 = Y written (4) even X0al. 0 ¯ have 0 0 to characterize inside X or outside X , existing −1 (X1 ), X1 = f −1 (X2 ), . . . , Xk ¯ interval methods will need to 0⊂bisect inside the 0 ⊃ box [u]k Mullier et = (2013). f −1 (X1 ),Equivalently, X1 = f −1 (X2we ), . could . . , Xk¯ have = Y written (4) even X0al. interval methods will need to bisect inside the box [u] ¯ to characterize inside X or outside X , existing k 0 0This is not (X2 ), . . .to , Xk¯solve = Y a chain (4) to amounts control the accuracy of the characterization. and the Xcomputation of1 = X0f −1 0 = f −1 (X1 ), X interval methods will need tocharacterization. bisect inside theThis boxis [u] ¯ k to solve a chain to control the accuracy of the not and the Xcomputation of1 = X0f amounts = f (X ), X (X ), . . . , X = Y (4) ¯ 0 1 2 k the problem interval methods will todimension. bisect inside theThis boxisthere [u] of settheinversion problems. Inamounts both cases, satisfactory due to itsneed large Note that to control the accuracy of the characterization. not to solve a chain and computation of X 0 of problems. In both cases, the problem satisfactory to its of large dimension. Note This that isthere to control thedue accuracy the characterization. not to solve a chain andset theinversion computation of X can be solved using a Set Inversion Jaulin and satisfactory exists numerical methods Mitchell (2003)Note to compute an 0Inamounts due to its large dimension. that there of set problems. both approach cases, the problem can be inversion solved using a Set Inversion approach Jaulin and exists numerical methods Mitchell (2003)Note to compute an satisfactory due to its large dimension. that there of set inversion problems. In both cases, the problem Walter (1993) based on interval analysis. Set inversion and exists outer approximation of theMitchell maximal backward reach set numerical methods (2003) to compute an can be solved using a Set Inversion approach Jaulin Walter (1993) based on interval analysis. Set inversion and outer approximation of the maximal backward reach set numerical methods (2003) to compute an can be (1993) solved based using on a Set Inversion approach Jaulin and and exists outer approximation of theMitchell maximal backward reach set Walter interval analysis. Set inversion outer approximation of the maximal backward reach set Walter (1993) based on interval analysis. Set inversion and Copyright © 2017 IFAC 9674 ∗ ∗ ∗ ∗
Copyright © 2017, 2017 IFAC 9674Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © IFAC (International Federation of Automatic Control) Copyright © 2017 IFAC 9674 Peer review under responsibility of International Federation of Automatic Copyright © 2017 IFAC 9674Control. 10.1016/j.ifacol.2017.08.1297
Proceedings of the 20th IFAC World Congress Benoit Desrochers et al. / IFAC PapersOnLine 50-1 (2017) 9266–9271 Toulouse, France, July 9-14, 2017
X⊃ 0 , or an inner approximation of the minimal backward reach set X⊂ 0 , but the main difficulty in inside the penumbra, i.e., to compute accurately, an inner approximation ⊂ of X⊃ 0 and an outer approximation X0 . The main contribution of this paper is to show that by rewriting the problem as a chain of set inversion problem, we will be able to characterize, with an arbitrary accuracy ¯ k the feasible set X0 , without any bisection of [u] . An ⊂ inner approximation of the penumbra X∂0 = X⊃ 0 \X0 , which corresponds to the uncertainty of the approximation, will also be computed. This was not possible before except for some specific cases such as when the system is linear Girard (2003). The paper is organized as follows. Section 2 formulates the set inversion problem in the case where it depends on a parameter vector p which is assumed to be inside a box [p]. Section 3 introduces the specific case where the function to be inverted is a link function. This will allow us to build, in Section 4, a large class of set inversion problems that can be solved efficiently, without any bisection inside the parameter space. Test-cases related to reachability analysis are presented in Section 5. Section 6 concludes the paper. 2. PROBLEM Consider the set inversion problem (8) X = f −1 (Y) n m where f : R → R , and Y is a subset of Rm . Solving a set inversion problem consists of characterizing the solution set X from inside and from outside. We consider the uncertain case where f depends on a parameter vector p ∈ Rq . The function should thus be written as f (x, p), but for simplicity of notation when we will compose functions to form a chain, we will often write fp (x). In such a case, the solution set X also depends on the value of p and is denoted by X(p). We have: (9) X⊂ ⊂ X(p) ⊂ X⊃ , where fp−1 (Y) = {x|∀p ∈ [p], fp (x) ∈ Y} X⊂ =
and will be called a thick set inversion problem. To characterize X with an arbitrary accuracy, i.e., to characterize the sets X⊂ and X⊃ , existing interval methods, will need to bisect with respect to the p-space, which makes them inefficient. Moreover, they are not able to prove that a box is inside the penumbra. It has been show recently Desrochers and Jaulin (2017) how it can proved that a box is inside the penumbra. It will be shown in Section 3 that, in the specific case where the function fp (x) is a link function (i.e., for all x, and for any box [p], the set f[p] (x) is a box), bisections on [p] can be avoided. We now illustrate how this can be done on a simple example. Example 1. Consider the thick set inversion problem −1 X = f[p] ([y]) where fp (x) =
p∈[p]
Denote by (P(Rn ), ⊂) the powerset of Rn equipped with the inclusion ⊂ as an order relation. The powerset P(Rn ) is a complete lattice with respect to ⊂. The meet operator corresponds to the intersection and the join to the union. The pair [X⊂ , X⊃ ] is an interval in P(Rn ) which is called a thick set and denoted by X . A thick set partitions Rn into three zones: the clear zone X⊂ , the penumbra X∂ = X⊃ \X⊂ and the dark zone Rn \X⊃ . A thick set X is a sub-lattice of (P(Rn ), ⊂), i.e., if A ∈ X ,B ∈ X , then A ∩ B ∈ X and A ∪ B ∈ X .
Notation. The set inversion problem we want to solve will be written as −1 X = f[p] (Y) (10)
p1 p3
p2 p4
x1 x2
.
(11)
− − + + p1− x1 + p2− x2 , p1+ x1 + p2+ x2 −p3 x1 +−p4 x 2 , p 3 x+1 + p4 x+2 p1 x1 + p2 x2 p 1 x 1 + p2 x 2 = , − + p− x + p x p+ 1 2 1 + p4 x 2 3 3 x 4 − + = f[p] (x), f[p] (x)
f (x, [p]) =
which corresponds to a box. If x ∈ [x] = [0, 1] × [2, 3], f − (x, [p]) belongs to the box 2 · [0, 1] + 2 · [2, 3] − ([x]) = f[p] 4 · [0, 1] − 6 · [2, 3] [0, 2] + [4, 6] = [0, 4] + [−18, −12] [4, 8] = [−18, −8]
fp−1 (Y) = {x|∃p ∈ [p], fp (x) ∈ Y} .
The set X∂ = X⊃ \X⊂ is called the penumbra Grimson and Lozano-P´erez (1985) and contains the unknown boundary of X.
Assume that [p] = [2, 3] × [2, 5] × [4, 5] × [−6, −1] and [y] = [5, 19] × [−7, 11]. Take a box [x] = [0, 1] × [2, 3] and let us show how is can be proved that it is inside the penumbra. Note that for simplicity, we have taken all intervals with a constant sign, but the method can be made much more general if we use modal interval analysis Goldsztejn et al. (2005) Vinas et al. (2006) or symbolic interval arithmetic Jaulin and Chabert (2010). Inside [x], we have
p∈[p]
X⊃ =
9267
(12)
and f + (x, [p]) belongs to the box 3 · [0, 1] + 5 · [2, 3] + ([x]) = f[p] 5 · [0, 1] − 1 · [2, 3] [0, 3] + [10, 15] = [0, 5] + [−3, −2] [10, 18] = . [−3, 3]
(13)
− + ([x]) ∩ [y] = ∅ and f[p] ([x]) ⊂ [y], we Since f[p] conclude that [x] is inside the penumbra. This is illustrated by Figure 1.
9675
Proceedings of the 20th IFAC World Congress 9268 Benoit Desrochers et al. / IFAC PapersOnLine 50-1 (2017) 9266–9271 Toulouse, France, July 9-14, 2017
0 -0.5
0
0.5
1
1.5
8
3.4
6
3.2
5
10
15
20
25
10
3.6
4 2
3 0 2.8 2.6 2.4 2.2 2
-2 -4 -6 -8 -10 -12
1.8 1.6 1.4
-14 -16 -18
+ Fig. 1. For [x] = [0, 1]×[2, 3] the set f[p] (blue dots) is inside − [y] whereas f[p] (red dots) is outside. We conclude that [x] is inside the penumbra
Fig. 2. Left: Classical interval methods accumulate on the yellow thick boundary (the penumbra). Right: the method we propose here allows a fast treatment of the penumbra. The frame box is [−5, 11]2 3. LINK FUNCTIONS
In the specific case where the function fp (x) is a link function (see Section 3 for more details), for all x, the set − + f[p] (x) is a box, denoted by f[p] (x), f[p] (x) , with lower
− + bound f[p] (x) ∈ Rm and upper bound f[p] (x) ∈ Rm . These two bounds can be described by an algorithm and two inclusion functions can be obtained for them using the rules of interval computation Kearfott (1996). This means that if x ∈ [x], we can obtain, under the form of an − + algorithm, two boxes f[p] ([x]) and f[p] ([x]) such that
− − + + f[p] ([x]) and f[p] ([x]). ([x]) ∈ f[p] ([x]) ∈ f[p] To characterize the thick set X of the thick set inversion −1 problem X = f[p] (Y), we start with a huge box (large enough to contain X⊃ ) and we cut it into several other boxes all stored inside a list L. For each box [x] of L, and in − + we compute the two boxes f[p] ([x]) and f[p] ([x]). From
these two boxes, we are able to test if [x] ⊂ X⊂ or if [x] ∩ X⊃ = ∅ or if [x] is included inside the penumbra. If nothing can be concluded, the box [x] is bisected and the two resulting boxes are stored inside the list. The treatment is applied to all boxes of the list until all boxes become too small (i.e., with a width smaller than a threshold ε). Note that, if we are able to conclude that a box [x] is inside the penumbra, no bisection would occur.
For Example 1, the method generates Figure 2 (right) which is an approximation of a thick set with the inner part (red), the outer part (blue) and the penumbra (orange). A classical interval method would not able to characterize the penumbra and yields Figure 2 (left). Since the resulting algorithm bisects in the x-sp ace, its complexity is exponential with respect to the dimension of x. Its application is thus limited to low dimensions. If do not bisect on the p-space (as we suggest in this paper using chains and link functions), the complexity with respect to p can be considered as polynomial.
This section considers the case of thick set-inversion problems in the case where the functions to be inverted have some good properties. These functions will be the links to build more complex functions to be inverted. Definition 2. Link function. If, for a box [p] ⊂ Rq , the set f[p] (x) = f (x, [p]) = {y ∈ Rp | ∃p ∈ [p] , y = fp (x)} is a box, then the function f is said to be a link. Links will be composed later in Section 4 in order to build a chain. If m = 1 then the function f (x, p) becomes scalar. Note that all f (x, p) that are scalar and continuous with respect to p, are links. Due to their specific box-shaped structure, link functions can be inverted with respect to x without bisections with respect to the uncertain parameter box [p]. Example 3. The function f (x, p) = 20e−x1 p − 8e−x2 p (14) is a scalar link function, due to the continuity of f with respect to p. Now, to get a computable expression for f (x, [p]), we need to compute the extremum of the function inside the interval [p] = [p− , p+ ]. Since ∂f (x, p¯) =0 ∂p ⇔ −20x1 · e−x1 p¯ + 8x2 · e−x2 p¯ = 0 2x2 (15) ⇔ e(x2 −x1 )p¯ = 5x 1 1 2x2 , ⇔ p¯ = ln x2 − x1 5x1 inside the interval [p] = [p− , p+ ], f (x, p), may have the extremum y¯ = f (x, p¯). As a result, if we define y − = f (x, p− ), y + = f (x, p+ ), we have − + min{y , y , y¯}, max{y − , y + , y¯} if p¯ ∈ [p] f (x, [p]) = if p¯ ∈ / [p] min{y − , y + }, max{y − , y + } Equivalently, − + (x), f[p] (x) (16) f (x, [p]) = f[p] where
9676
− (x) f[p]
=
min{y − , y + } if p¯ ∈ / [p] min{y − , y + , y¯} if p¯ ∈ [p]
(17)
Proceedings of the 20th IFAC World Congress Benoit Desrochers et al. / IFAC PapersOnLine 50-1 (2017) 9266–9271 Toulouse, France, July 9-14, 2017
and + (x) f[p]
=
/ [p] max{y − , y + } if p¯ ∈ max{y − , y + , y¯} if p¯ ∈ [p]
9269
(18)
− + (x) and f[p] (x) belong to R. If x ∈ [x], Note that both f[p] we have − + p] ∩ [p] = ∅ min{ − − −y , +y } if [¯ ([x]) = (x) ∈ f[p] f[p] y ]} otherwise min{ y , y , [¯
and
− + p] ∩ [p] = ∅ max{ −y , +y } if [¯ y ]} otherwise Fig. 3. Left: measurement boxes [t ] × [y ]. Right: Approxmax{ y , y , [¯ i i where imation of the thick set Q in the parameter space 1 2 [x2 ] ln [¯ p] = Example 7. Consider the problem of estimating the pa[x2 ] − [x1 ] 5 [x1 ] y] = f ([x] , [¯ p] ∩ [p]) (19) rameters q1 and q2 of the model [¯ − − y(q, t) = 20e−q1 t − 8e−q2 t . f ([x] , p ) y + = y = f ([x] , p+ ). We assume that 10 measurements yi have been collected at Definition 4. Vector components. The family of vectors time ti . The uncertainties on the pair (ti , yi ) is represented {p1 , p2 , . . . , p } correspond to vector components of the by intervals as represented by the table below and Figure vector p if each pi is a subvector of p and the indexes are 3 (left). We are interested by the thick set all disjoint. For instance if p = (p1 , . . . , p9 ), then {p1 , p2 }, Q = q ∈ R2 |∀i, y(q, ti ) ∈ [yi ] . where p1 = (p2 , p9 ) and p2 = (p8 , p1 , p3 ), correspond to If we set vector components of p. Note that the two index sets {2, 9} y(q, [t1 ]) and {8, 1, 3} are disjoint. .. f[t] (q) = . The fact that index sets are disjoint implies in independeny(q, [t10 ]) cies between parameters. The conservatism related to the [y] = [y10 ] × · · · × [y10 ] dependency effect, well known in interval analysis Moore [t] = [t1 ] × · · · × [t10 ] (1979), is thus avoided. In practice, this will allow us to then, we have avoid bisection in the {p1 , p2 , . . . , p }-space. −1 Q = f[t] ([y]), Proposition 5. Consider scalar link functions, fi (x, pi ) , where vectors {p1 , p2 , . . . , p } are vector components of which is thick set inversion problem with the link function the vector p. The function f[t] (q) to be inverted. The inversion yields the approxima tion of Figure 3 (right). f1 (x, p1 ) . .. (20) f (x, p) = i [ti ] [yi ] 1 [0.25, 1.25] [2.7,12.1] f (x, p ) 2 [1, 2] [1.04,7.14] is a link. 3 [1.75, 2.75] [-0.13,3.61] 4 [2.5, 3.5] [-0.95,1.15] Proof. Consider the box [p] and its box components 5 [5.5, 6.5] [-4.85,-0.29] [p1 ] , . . . , [p ]. The set f (x, [p]) is the Cartesian product 6 [8.5, 9.5] [-5.06,-0.36] of intervals: 7 [12.5, 13.5] [-4.1,-0.04] f (x, [p]) = f1 (x, [p1 ]) × · · · × f (x, [p ]) (21) 8 [16.5, 17.5] [-3.16,0.3] 9 [20.5, 21.5] [-2.5,0.51] which corresponds to a box. 10 [24.5, 25.5] [-2,0.6] Example 6. The function p1 x1 x2 + p1 p2 sin(x2 ) 4. CHAIN f (x, p) = (22) p3 x1 + p3 p4 cos(x1 )x2 + f[p] (x)
+ ∈ f[p] ([x]) =
is a link. Now, if we replace p3 by p2 , a dependency occurs between the two components of f and the function is not a link anymore.
Consequence. If fp (x) is a link, for all x, the set f[p] (x) is a box. If x ∈ [x], the lower and upper bounds of the − box f[p] (x) is included inside the boxes f[p] ([x]) and + f[p] ([x]). As a consequence, the methodology explained in Section 2 can be applied to find an inner and outer −1 approximations of the thick set X = f[p] (Y), with an approximation of the penumbra.
Definition 8. Chain. The function ϕ(x, p) = ϕp (x) is said to be a chain if it can be written as (23) ϕp (x) = fp ◦ · · · ◦ fp22 ◦ fp11 (x)
where {p1 , p2 , . . . , p } are vector components of p and the fpkk are links.
Note that the parameter vectors p1 , p2 , . . . , p are independent. This implies that the pessimism due to the dependency effect will not take place during the resolution. Moreover, the fact that fpkk are links implies that pessimism due to the wrapping effect (of a set by a box) will not exist.
9677
Proceedings of the 20th IFAC World Congress 9270 Benoit Desrochers et al. / IFAC PapersOnLine 50-1 (2017) 9266–9271 Toulouse, France, July 9-14, 2017
Example 9. The function
sin(p3 x21 ) + p4 x ϕp (x) = p5 (p1 e 2 + p1 p2 x1 x2 ) p3 x21
(24)
is not a link. Indeed, p3 occurs on both components of ϕp . Now, if we define p 1 e x 2 + p 1 p2 x 1 x 2 1 fp1 (x) = (25) p3 x21 fp22 (z)
=
sin(z2 ) + p4 p 5 z1 z2
(26)
with p1 = (p1 , p2 , p3 ) and p2 = (p4 , p5 ) then, we have ϕp (x) = fp22 ◦ fp11 (x).
(27)
Since both fp22 , fp11 are links, ϕp (x) is a chain. Consequence. Solving a thick set inversion problem X = ϕ−1 [p] (Y) in the case where ϕ is a chain can be done by solving several thick set inversion problems, without any bisection inside the box [p]. Rewriting a chain as a composition of link functions can probably be done by symbolic methods, but we do not know any approach or algorithm to perform this task. Of course, we should not allow the dimension of the intermediate spaces to be arbitrary large compare to the dimension of x if we want to be efficient. The following section provides some test-cases of the thick inversion of chains, applied to reachability problems.
5. TEST-CASES
Fig. 4. Approximation of the sets X(0), X(1), X(2), X(3) Test-case 2. Consider the nonlinear system x (k) + x2 (k) · u (k) 1 1 2 x1 (k + 1) = 1 (31) x2 (k + 1) · x1 (k) · x2 (k) + u2 (k) 2 with u1 ∈ [1, 2] and u2 ∈ [−2, −1]. The set Y is assumed to be a centred disk which has to be reached at time k¯ = 3. Solving the chain of thick set inversion problem generates Figure 5. The frame boxes are [−12, 10] × [−12, 10] for all X (i).
Test-case 1. Consider the linear system x(k + 1) = A · x(k) where A(k) ∈
[2.5, 3] [2, 3] [4, 4.5] [−3, −2]
(28)
(29)
plays the role of the parameter vector p. Assume that we want to reach the target Y = [4, 20] × [−8, 12] (30) ¯ at time k = 3. Solving the chain of thick set inversion problems yields the approximation of the thick sets X (0), X (1), X (2), X (3) as depicted on Figure 4. Note that X (3) corresponds to the box Y. We first compute the thick set X (2) = f −1 [A] ( X (3)). Note that X (3) is here a thin set (i.e., a classical set) whereas X (2) is thick (i.e., with a penumbra). Then we compute X (1) = f −1 [A] ( X (2)) and finally, we compute X (0) = −1 f [A] ( X (1)). At each step, the penumbra inflates, due to the accumulation of uncertainties. The frame boxes are [−0.6, 0.8] × [−0.4, 1] for X (0), [−1.5, 3]×[−2.5, 2.2] for X (1),[−4, 8]×[−3, 10] for X (2), and [−5, 27] × [−15, 20] for X (3).
Fig. 5. Approximation of the sets X(0), X(1), X(2), X(3) 6. CONCLUSION In this paper, we have presented a new interval based approach to compute the set X0 of initial states that will
9678
Proceedings of the 20th IFAC World Congress Benoit Desrochers et al. / IFAC PapersOnLine 50-1 (2017) 9266–9271 Toulouse, France, July 9-14, 2017
reach a target Y in a finite time. The dynamic system that is considered is discrete time and uncertain (the uncertainties are represented by intervals). The principle is to transform the problem into a chain of set inversion problems composed with links, i.e., a function with nice properties for set inversion. Among the problems that remain to be solved, we need (1) to find symbolic methods for the decomposition of a chain into links with low dimensions; (2) to extend the approach to continuous systems described by differential equations and to study how the penumbra can be characterized from inside in such a case; and (3) to solve the chain of set inversion problems in parallel and not sequentially as done here. REFERENCES Aubin, J.P. (2001). Viability Kernels and Capture Basins of Sets Under Differential Inclusions. SIAM Journal on Control and Optimization, 40(3), 853–881. doi: 10.1137/S036301290036968X. Aubin, J. and Frankowska, H. (1990). Set-Valued Analysis. Birkh¨ auser, Boston. Bayen, A.M., Mitchell, I.M., Osihi, M.K., and Tomlin, C.J. (2007). Aircraft Autolander Safety Analysis Through Optimal Control-Based Reach Set Computation. Journal of Guidance, Control, and Dynamics, 30(1), 68–77. doi:10.2514/1.21562. Blanchini, F. and Miani, S. (2007). Set-Theoretic Methods in Control. Springer Science & Business Media. Bouissou, O., Chapoutot, A., Djaballah, A., and Kieffer, M. (2014). Computation of parametric barrier functions for dynamical systems using interval analysis. In IEEE CDC. Los Angeles, United States. Desilles, A., Zidani, H., and Cruck, E. (2012). Collision analysis for an UAV. In AIAA Guidance, Navigation, and Control Conference, 13–16. American Institute of Aeronautics and Astronautics. Desrochers, B. and Jaulin, L. (2017). Computing a guaranteed approximation the zone explored by a robot. IEEE Transaction on Automatic Control, 62(1), 425– 430. dit Sandretto, J.A., Trombettoni, G., Daney, D., and Chabert, G. (2014). Certified calibration of a cabledriven robot using interval contractor programming. In F. Thomas and A.P. Gracia (eds.), Computational Kinematics, Mechanisms and Machine Science, Spinger. Drevelle, V. and Bonnifait, P. (2012). igps: Global positioning in urban canyons with road surface maps. IEEE Intelligent Transportation Systems Magazine, 4(3), 6– 18. Girard, A. (2003). Computation and stability analysis of limit cycles in piecewise linear hybrid systems. In 1st IFAC Conference on Analysis and Design of Hybrid Systems, 181–186. Goldsztejn, A., Daney, D., Rueher, M., and Taillibert, P. (2005). Modal intervals revisited: a mean-value extension to generalized intervals. In Proceedings of the Quantification in Constraint Programming, Barcelona, Spain, 2005 (QCP 2005). Grimson, W.E.L. and Lozano-P´erez, T. (1985). Recognition and localization of overlapping parts in two and three dimensions. In IEEE International Conference on Robotics and Automation, 61–66.
9271
Jauberthie, C., Verdi`ere, N., and Trav´e-Massuy`es, L. (2011). Set-membership identifiability: definitions and analysis. In In Proceedings of the 18th IFAC World Congress, Milan, Italie, volume 00, 12024–12029. aaa. Jaulin, L. and Chabert, G. (2010). Resolution of nonlinear interval problems using symbolic interval arithmetic. Engineering Applications of Artificial Intelligence, 23(6), 1035–1049. Jaulin, L. and Walter, E. (1993). Guaranteed nonlinear parameter estimation via interval computations. Interval Computation, 61–75. Jaulin, L. and Walter, E. (1994). Modeling Techniques for Uncertain Systems, chapter Set Inversion, with Application to Guaranteed Nonlinear Estimation and Robust Control, 3–20. Birkh¨auser. Kearfott, R.B. (1996). Rigorous Global Search: Continuous Problems. Kluwer, Dordrecht, the Netherlands. Konecny, M., Taha, W., Duracz, J., Duracz, A., and Ames, A. (2013). Enclosing the behavior of a hybrid system up to and beyond a zeno point. In Cyber-Physical Systems, Networks, and Applications (CPSNA). Mitchell, I.M. (2003). Application of level set methods to control and reachability problems in continuous and hybrid systems. Moore, R.E. (1979). Methods and Applications of Interval Analysis. SIAM, Philadelphia, PA. Mullier, O., Goubault, E., Kieffer, M., and Putot, S. (2013). General inner approximation of vector-valued functions. Reliable Computing, 18(0), 117–143. Raissi, T., Efimov, D., and Zolghadri, A. (2012). Interval state estimation for a class of nonlinear systems. IEEE Transactions on Automatic Control, 57(1), 260–265. Raissi, T., Ramdani, N., and Candau, Y. (2004). Set membership state and parameter estimation for systems described by nonlinear differential equations. Automatica, 40, 1771–1777. Ramdani, N. and Nedialkov, N. (2011). Computing Reachable Sets for Uncertain Nonlinear Hybrid Systems using Interval Constraint Propagation Techniques. Nonlinear Analysis: Hybrid Systems, 5(2), 149–162. Taha, W. and Duracz, A. (2015). Acumen: An opensource testbed for cyber-physical systems research. In CYCLONE’15. Vinas, P.H., Sainz, M.A., Vehi, J., and Jaulin, L. (2006). Quantified set inversion algorithm with applications to control. Reliable computing, 11(5), 369–382.
9679