Robust-to-Dynamics Optimization

Robust-to-Dynamics Optimization ∗ Amir Ali Ahmadi † and Okta y G¨ unl¨ uk ‡ Abstract. A r obust-to-dynamics optimization (RDO) pr oblem is an optimization problem sp eciﬁed by tw o pieces of input: (i) a mathematical program (an ob jective function f : R n → R and a feasible set Ω ⊆ R n ), and (ii) a dynamical system (a map g : R n → R n ). Its goal is to minimize f ov er the set S ⊆ Ω of initial conditions that forever remain in Ω under g . The focus of this pap er is on the case where the mathematical program is a linear program and the dynamical system is either a kno wn linear map, or an uncertain linear map that can change ov er time. In b oth cases, we study a conv erging sequence of p olyhedral outer approximations and (lifted) spectrahedral inner appro ximations to S . Our inner appro ximations are optimized with resp ect to the ob jective function f and their semideﬁnite c haracterization—whic h has a semideﬁnite constrain t of ﬁxed size— is obtained by applying p olar duality to conv ex sets that are inv arian t under (multiple) linear maps. W e characterize three barriers that can stop conv ergence of the outer appro ximations to S from b eing ﬁnite. W e pro ve that once these barriers are remo ved, our inner and outer appro ximating pro cedures ﬁnd an optimal solution and a certiﬁcate of optimalit y for the RDO problem in a ﬁnite num b er of steps. Moreov er, in the case where the dynamics are linear, w e show that this phenomenon o ccurs in a num b er of steps that can b e computed in time p olynomial in the bit size of the input data. Our analysis also leads to a p olynomial-time algorithm for RDO instances where the sp ectral radius of the linear map is b ounded abov e by any constant less than one. Finally , in our concluding section, w e prop ose a broader research agenda for studying optimization pr oblems with dynamic al systems c onstr aints , of which RDO is a sp ecial case. Key wo rds. Optimization in dynamical systems, semi-inﬁnite linear programs, joint sp ectral radius, semideﬁ- nite programming-based approximations. AMS subject classiﬁcations. 90C34, 90C25, 90C22, 93D05, 93B40. 1. Intro duction. In man y real-w orld situations, a decision mak er is called upon to mak e a decision that optimizes an ob jectiv e function and satisﬁes a current set of constraints. These constrain ts ho w ev er can c hange o v er time b ecause of the inﬂuence of an external dynamical system, rendering the original decision infeasible. F or example, a v accine that meets a ce rtain lev el of eﬃcacy against a given virus can fail to do so under future m utations of the virus. The goal is then to make the b est decision among all p ossible decisions which remain feasible as the constrain ts v ary o v er time. By changing the p oin t of reference, suc h a scenario can equiv alently b e thought of as a situation where a constrain t set is ﬁxed, but the original decision is mov ed around under the inﬂuence of a dynamical system. The goal is then to mak e the b est decision among all decisions that remain in the (ﬁxed) constraint set under the inﬂuence of the dynamical system. ∗ This wo rk w as partially funded b y the DARP A Y oung Facult y Awa rd, the Y oung Investigator Awa rd of the AF OSR, the CAREER Awa rd of the NSF, the Go ogle Facult y Awa rd, and the Sloan Fello wship. The authors w ould also like to thank the Mathematical Institute at Oxford University for hosting them while part of this resea rch was b eing conducted. † ORFE, Princeton University , Sherrerd Hall, Princeton, NJ 08540 ( aaa@princeton.edu ) ‡ Co rnell University , Scho ol of ORIE, Ithaca, NY 14853 ( ong5@cornell.edu ). The author was partially supp orted b y ONR grant N00014-21-1-2575 1 2 A. A. AHMADI AND O. G ¨ UNL ¨ UK In this pap er, w e study a mathematical abstraction of problems of this nature. W e will refer to them as “r obust-to-dynamics optimization ” (RDO) problems. The name alludes to the fact that the solution to these problems needs to b e robust to external dynamics, in the sense that it should remain feasible at all times as it is mo v ed by a dynamical system. RDO problems hav e a very natural mathematical formulation which relies on tw o pieces of input: 1. an optimization problem: (1.1) min x { f ( x ) : x ∈ Ω } , 2. a dynamical system: (1.2) ( x k +1 = g ( x k ) in discrete time (fo cus of this pap er), or ˙ x = g ( x ) in contin uous time . Here, we hav e x ∈ R n , Ω ⊆ R n , f : R n → R , g : R n → R n ; x k denotes the state at time step k , and ˙ x is the deriv ative of x with resp ect to time. RDO is then the following optimization problem: 1 (1.3) min x 0 { f ( x 0 ) : x k ∈ Ω for k = 0 , 1 , 2 , . . . , u.t.d. x k +1 = g ( x k ) } in discrete time, or min x 0 { f ( x 0 ) : x ( t ; x 0 ) ∈ Ω , ∀ t ≥ 0 , u.t.d. ˙ x = g ( x ) } in con tin uous time, where x ( t ; x 0 ) denotes the solution of the diﬀerential equation ˙ x = g ( x ) at time t , starting at the initial condition x 0 ∈ R n . By the u.t.d. notation, w e imply that x k (resp. x ( t ; x 0 )) must satisfy the equation of our dynamics x k +1 = g ( x k ) (resp. ˙ x = g ( x )) for k = 0 , 1 , 2 , . . . (resp. t ≥ 0). In w ords, we are optimizing an ob jective function f ov er the set of initial conditions that never lea ve the set Ω under the dynamics gov erned b y g . RDO problems can naturally b e categorized dep ending on the type of optimization prob- lem considered in ( 1.1 ) and the t yp e of dynamics considered in ( 1.2 ). A list of some p ossible com binations is giv en in the table b elow. Optimization Problem ( f , Ω) Dynamical System ( g ) Linear program Linear Con v ex quadratic program Nonlinear (e.g., p olynomial) Second order cone program Uncertain Semideﬁnite program Time-v arying In teger program Uncertain and time-v arying . . . . . . T able 1: Any combination of an en try from the ﬁrst column and an en try from the second column of this table leads to an RDO problem. 1 W e write “u.t.d.” as an abbreviation for “under the dynamics”. ROBUST-TO-D YNAMICS OPTIMIZA TION 3 This framew ork can b e even further generalized b y considering optimization problems with “dynamical system constraints” (see Section 4 ). In this pap er, w e fo cus on RDO problems where the mathematical program is a linear program and the dynamical system is discrete and either a known linear map, or an uncertain linear map that can change ov er time. The second case subsumes some uncertain nonlinear dynamics that are not time-v arying as well. As is made eviden t ab o v e, RDO is at the juncture b etw een optimization and dynamical systems, and as suc h, has some commonalities with literature in b oth areas. On the opti- mization side, RDO comes closest to the area of robust optimization and on the dynamical systems side, it has connections to the theory of inv arian t sets. W e describ e the links b etw een these areas in more detail b elo w. In its most common form, robust optimization (R O) [ 9 , 12 , 39 , 13 , 10 ] deals with problems of the t yp e (1.4) min x { f ( x ) : u i ( x ) ≤ 0 ∀ u i ∈ U i , i = 1 , . . . , m } , where U i is a prescrib ed uncertaint y set for (the parameters of ) the constrain t function u i : R n → R . Lik e RDO, RO problems ha v e to contend with uncertain t y in the constraints, though unlike RDO, this uncertain t y is not explicitly linked to a dynamical system. As an area, RO is well-studied from a computational complexity standp oin t. By now, w e almost fully understand when a robust optimization problem inv olving a particular mathematical program (e.g., a linear or a conv ex quadratic program) and a particular type of uncertaint y set U i (e.g., p olyhedral, ellipsoidal, etc.) is p olynomial-time solv able or NP-hard; see [ 12 ] and [ 9 ] and references therein. W e note that similar to robust optimization problems, the RDO framew ork is a sp ecial case of semi-inﬁnite optimization [ 36 , 23 , 43 ]. In the case of RDO problems with discrete deterministic dynamics, one is dealing with a countable semi-inﬁnite optimization problem. These constraints, ho w ev er, hav e a very sp ecial structure arising from the underlying dynamical system, and we exploit this structure hea vily in this pap er. On the dynamics side, inv ariant set theory [ 16 , 32 , 15 , 25 , 38 , 33 , 7 ] concerns itself with the study of sets that are inv ariant under the action of dynamical systems. It also considers the problem of designing con trollers that inﬂuence a dynamical system so as to mak e a desired set inv ariant. This problem has applications to mo del predictiv e control [ 31 , 24 , 35 ], among other subﬁelds in control theory . The literature on inv ariance in control b y and large do es not consider constraints and studies the existence and the structure of inv arian t sets for diﬀeren t t yp es of dynamical systems. The subset of the literature that do es incorp orate constrain ts t ypically aims at characterizing the maximal (with resp ect to inclusion) inv ariant set within a giv en constrain t set. These maximal inv ariant sets are often complicated to describ e ev en for simple dynamical systems and hence are appro ximated. While both inner and outer appro ximations are interesting, inner appro ximations are more relev ant to applications as they provide initial conditions that remain within the constrain t set for all time. T o the b est of our kno wledge, the appro ximations a v ailable in the literature do not take in to consideration an ob jectiv e function as is done in this pap er (cf. Section 2.2.2 and Section 3.2.2 ). While problems related to RDO hav e b een studied in the control communit y , we b elieve that the framew ork in whic h w e presen t these problems—as a merging of a mathematical program and a dynamical system (cf. T able 1 )—and the questions that w e study , and hence 4 A. A. AHMADI AND O. G ¨ UNL ¨ UK our results, are diﬀeren t. Our hop e is that this framew ork will make problems of this type more palatable to the optimization communit y . W e also b elieve that the framew ork provides the right setup for a systematic algorithmic study of RDO, as has b een done so successfully for RO. Indeed, our ov erarching goal is to pro vide an understanding of the computational complexit y of each type of RDO problem that arises from the tw o columns of T able 1 . This can b e either in the form of negative results (e.g., NP-hardness or undecidabilit y pro ofs), or p ositiv e ones (algorithms with guarantees). The current pap er is a step in this direction. V ery recently , the work presen ted in this pap er has b een applied to the problem of safely learning dynamical systems from tra jectory data [ 1 ]. This problem has applications to rein- forcemen t learning and rob otics and requires the learning pro cess to resp ect certain “safety constrain ts”. These constraints restrict the set of allo w able queries (i.e., starting p oin ts) to the unknown dynamical system to those whose future tra jectories are guaranteed to remain in a given (safe) set; see [ 1 ] for details. W e also b eliev e there are p otential applications to constrained nonlinear optimization. Indeed, one can think of any algorithm for minimizing a function f : R n → R (e.g., gradien t descent) as a dynamical system. If the function is to b e minimized o v er a feasible set Ω ⊆ R n , then characterizing the set of initial conditions that remain feasible under the iterations of the algorithm is an RDO problem. In the case where f is a quadratic form, Ω is p olyhedral, and the minimization algorithm is gradient descent (p ossibly with time-v arying step sizes), the resulting RDO problem falls within those studied in this pap er. 1.1. Outline and contributions of the pap er. In Section 2 , we study robust to linear dynamics linear programs (R-LD-LPs), which are RDO problems where the optimization problem is a linear program and the dynamics are linear. W e show that the feasible set S of an R-LD-LP is not alwa ys p olyhedral and that testing membership of a p oint to S is NP-hard (Theorem 2.1 ). In Section 2.1 , we study a sequence of natural outer approximations S r to S that get tighter and tigh ter as r gro ws. W e giv e a p olynomial-time chec k able criterion for testing whether S r = S for a given nonnegative integer r (Lemma 2.3 ). W e then c haracterize three conditions that ma y stop con v ergence of S r to S from happ ening in a ﬁnite num ber of steps (Propositions 2.4 , 2.5 , 2.6 , and Theorem 2.1 , part(ii)). These conditions hav e previously app eared as assumptions in the literature on inv ariant sets (see, e.g., [ 7 , Sect. 2.2] and references therein), but to the b est of our kno wledge, there are no formal arguments that sho w why all three of the conditions are needed to guarantee ﬁnite conv ergence. Once this is clariﬁed, the main theorem of Section 2.1 (Theorem 2.7 ) shows that under these three conditions, con v ergence of S r to S is not only ﬁnite but takes a n um b er of steps that can b e computed in time p olynomial in the size of the data. Our pro of also sho ws that all instances of R-LD-LP for which the sp ectral radius of the matrix deﬁning the linear dynamics is upp er b ounded by a constan t less than one can b e solved in p olynomial time. In Section 2.2 , we study inner appro ximations of S which hav e the adv an tage (compara- tiv ely to outer appro ximations) of pro viding feasible solutions to an R-LD-LP . W e ﬁrst give a general construction that starts with an y full-dimensional and compact inv ariant set for the dynamical system and pro duces a sequence of nested inner approximations to S that conv erge to S ﬁnitely (Lemma 2.8 , Lemma 2.9 , and Corollary 2.10 ). Using this pro cedure, a ﬁnitely- con v ergen t sequence of upp er b ounds on the optimal v alue of an R-LD-LP can b e computed ROBUST-TO-D YNAMICS OPTIMIZA TION 5 b y solving conv ex quadratic programs (cf. Section 2.2.1 ). In Section 2.2.2 , w e formulate a sequence of semideﬁnite programs whic h ﬁnd inner appro ximations to S that are optimally aligned with the ob jective function of the R-LD-LP (Theorem 2.11 , part(i)). W e sho w that the solutions to our semideﬁnite programs coincide with the optimal solutions to our R-LD-LP after a n um b er of steps that can b e computed in p olynomial time (Theorem 2.11 , part (ii)). In Section 3 , we study robust to unc ertain and time-varying linear dynamics linear pro- grams (R-UTVLD-LPs), which are RDO problems where the optimization problem is a linear program and the dynamical system is linear, but it is uncertain and time-v arying. This problem class also captures certain nonlinear time-inv ariant dynamics in the presence of un- certain t y . The algorithmic questions get considerably more in v olv ed here; e.g., ev en testing mem b ership of a given p oin t to the feasible set of an R-UTVLD-LPs is undecidable. The goal of the section is to generalize (to the extent p ossible) the results of Section 2 , b oth on outer and inner appro ximations of the feasible set. T o do this, we replace the notion of the sp ectral radius with that of the joint sp ectral radius (see Deﬁnition 3.1 ), and ellipsoidal in v ariant sets with inv ariant sets that are a ﬁnite in tersection of ellipsoids (cf. Theorem 3.4 ). W e prov e ﬁnite con v ergence of b oth inner and outer appro ximations and giv e a semideﬁnite form ula- tion of inner approximations that are aligned with the ob jectiv e function (cf. Lemma 3.3 , Theorems 3.2 and 3.6 ). Finally , in Section 3.3 , we provide a p olynomial-size extended LP form ulation of R-UTVLD-LP problems where the dynamical system is describ ed by the set of p ermutation matrices. W e conclude in Section 4 b y noting that RDO can b e inscrib ed in a broader class of prob- lems, those that w e refer to as “optimization problems with dynamical system constraints”. W e b eliev e that problems of this type form a rich class of optimization problems and hop e that our pap er will instigate further research in this area. 2. Robust to linear dynamics linea r programming. W e deﬁne a r obust to line ar dynamics line ar pr o gr am (R-LD-LP) to b e an optimization problem of the form 2 (2.1) min x 0 ∈ R n { c T x 0 : x k ∈ P for k = 0 , 1 , 2 , . . . , u.t.d. x k +1 = Gx k } , where P : = { x ∈ R n | Ax ≤ b } is a p olyhedron and G ∈ R n × n is a giv en matrix. One can equiv alently formulate an R-LD-LP as a linear program of a particular structure with an inﬁnite num b er of constraints: (2.2) min x ∈ R n { c T x : G k x ∈ P for k = 0 , 1 , 2 , . . . } . The input to an R-LD-LP is fully deﬁned 3 b y c ∈ R n , A ∈ R m × n , b ∈ R m , G ∈ R n × n . Note that a linear program can b e thought of as a sp ecial case of an R-LD-LP where the matrix G is the identit y matrix. Problem ( 2.1 ), or equiv alen tly problem ( 2.2 ), has a simple geometric in terpretation: we are in terested in optimizing a linear function not o v er the entire polyhedron 2 The dep endence of the ob jectiv e function on just the initial condition is without loss of generality . Indeed, if the ob jectiv e instead read P N k =0 ˆ c T k x k , w e could let c T = P N k =0 ˆ c T k G k in ( 2.1 ) to hav e an equiv alen t problem. 3 Whenev er we study complexity questions around an R-LD-LP , we use the standard T uring mo del of computation (see, e.g., [ 49 ]) and consider instances where the en tries of c, A, b, G are rational num b ers and hence the input can b e represented with a ﬁnite n umber of bits. 6 A. A. AHMADI AND O. G ¨ UNL ¨ UK P , but o v er a subset of it that do es not leav e the p olyhedron under the application of G , G 2 , G 3 , etc. Hence, the feasible set of an R-LD-LP is by deﬁnition the following set (2.3) S : = ∞ \ k =0 { x ∈ R n | AG k x ≤ b } = ∞ \ k =0 m \ i =1 { x ∈ R n | a T i G k x ≤ b i } where a T i and b i denote the i th row of A and b , resp ectiv ely . When written in this form, it is clear that R-LD-LP is a sp ecial case of a linear semi-iniﬁte program (see, e.g., [ 36 , Sect. 4]). In this section, w e exploit the sp ecial dynamics-based structure of R-LD-LP to characterize its complexity and determine sub classes of it that can b e solv ed in p olynomial time. W e start with a theorem on the basic prop erties of the set S in ( 2.3 ). W e recall that a set T ⊆ R n is said to b e invariant under a dynamical system, if all tra jectories of the dynamical system that start in T remain in T forever. W e will simply sa y that a set is inv ariant if the underlying dynamical system is clear from the con text. Theo rem 2.1. The set S in ( 2.3 ) has the fol lowing pr op erties: (i) It is close d, c onvex, and invariant. (ii) It is not always p olyhe dr al (even when A, b, G have r ational entries). (iii) T esting memb ership of a given p oint to S is NP-har d. Pr o of. W e prov e the three statemen ts separately . (i) Con v exit y and closedness are a consequence of the fact that p olyhedra are con v ex and closed, and that (inﬁnite) in tersections of conv ex (resp. closed) sets are conv ex (resp. closed). In v ariance is a trivial implication of the deﬁnition: if x ∈ S , then Gx ∈ S . (ii) Clearly the set S can b e p olyhedral (consider, e.g., the case where G = I , i.e., the iden tit y matrix). W e give an example where it is not. Consider an R-LD-LP in R 2 where P = [ − 1 , 1] 2 and (2.4) G =  4 / 5 3 / 5 − 3 / 5 4 / 5  =  cos( θ ) sin( θ ) − sin( θ ) cos( θ )  , with θ = arcsin(3 / 5). This is the smallest angle (in radians) of the right triangle with sides 3, 4, and 5. F rom Niven’s theorem [ 40 , Corollary 3.12], if ¯ θ ∈ [0 , π / 2] and sin( ¯ θ ) is rational, then ¯ θ /π is not rational unless ¯ θ ∈ { 0 , π / 6 , π / 2 } . Consequently , the num b er θ /π is irrational, whic h means that θ and π are rationally indep enden t. Also note that G is a rotation matrix that rotates points in the xy -Cartesian plane clo ckwise by θ . In other w ords, using p olar co ordinates with the conv ention x k = ( r k , ϕ k ), the feasible region of this R-LD-LP consists of p oin ts x 0 = ( r 0 , ϕ 0 ) such that ( r 0 , ϕ 0 + k θ ) ∈ P for all integers k ≥ 0. Notice that the closed disk D cen tered at the origin with radius 1 is contained in P and consequen tly ( r 0 , ϕ 0 + kθ ) ∈ P for all k pro vided that r 0 ≤ 1. On the other hand, consider a p oint x 0 = ( r 0 , ϕ 0 ) ∈ P such that r 0 > 1 . Let C b e the circle centered at the origin with radius r 0 . Clearly , x k ∈ C , ∀ k . F urthermore, note that for any ﬁxed r 0 > 1, there exist scalars β 1 , β 2 , with 0 ≤ β 1 < β 2 < 2 π , suc h that none of the p oints in C on the arc b et w een ( r 0 , β 1 ) and ( r 0 , β 2 ) b elong to P . Ho w ev er, as θ and 2 π are rationally indep endent, w e must ha v e ϕ 0 + k θ ∈ ( β 1 , β 2 ) for some integer k ≥ 1 (see, e.g., [ 21 , Chapter 3, Theorem 1]). Consequen tly , x 0 is feasible if and only if r 0 ≤ 1, i.e., S = D . ROBUST-TO-D YNAMICS OPTIMIZA TION 7 If w e let S r : = T r k =0 { x ∈ R 2 | G k x ∈ [ − 1 , 1] 2 } , Figure 1 depicts that as r increases, more and more p oin ts lea v e the p olytope P , until nothing but the unit disk D is left. S 0 S 2 S 4 S 16 Figure 1: The construction in the pro of of Theorem 2.1 , part (ii). (iii) W e no w pro v e that the follo wing decision problem is NP-hard ev en when m = 1: Giv en z ∈ Q n , A ∈ Q m × n , b ∈ Q m , G ∈ Q n × n , test if z / ∈ S . W e show this via a p olynomial- time reduction from the following decision problem, which is known to b e NP-hard (see [ 18 , Corollary 1.2]): Giv en a directed graph Γ on n no des, test if there exists an integer ˆ k ≥ 1 such that Γ has no directed paths 4 of length ˆ k starting from no de 1 and ending at no de n . Let G ∈ Q n × n b e the adjacency matrix of the graph Γ, i.e., a matrix with its ( i, j )-th en try equal to 1 if there is an edge from no de i to no de j in Γ and equal to 0 otherwise. Let z = (0 , . . . , 0 , 1) T ∈ Q n , A = ( − 1 , 0 , . . . , 0) G ∈ Q 1 × n , and b = − 1 2 . Let S b e as in ( 2.3 ). W e claim that z / ∈ S if and only if for some integer ˆ k ≥ 1, Γ has no directed paths of length ˆ k from no de 1 to no de n . T o argue this, recall that the ( i, j )-th entry of G k is greater than or equal to one if and only if there exists a directed path of length k in Γ from no de i to no de j . Supp ose ﬁrst that for all integers k ≥ 1, Γ has a directed path of length k from no de 1 to no de n. This implies that (1 , . . . , 0) G k (0 , . . . , 1) T ≥ 1 , for k = 1 , 2 , 3 , . . . , which in turn implies that AG k − 1 z ≤ − 1 for k = 1 , 2 , . . . . As b = − 1 2 , we ha v e AG k z ≤ b , for k = 0 , 1 , . . . , and hence z ∈ S . Supp ose next that for some integer ˆ k ≥ 1, Γ has no directed paths of length ˆ k from no de 1 to no de n . This implies that (1 , . . . , 0) G ˆ k (0 , . . . , 1) T = 0 and hence w e ha v e AG ˆ k − 1 z = 0 > b . As z / ∈ S ˆ k − 1 , z / ∈ S . P art (iii) of Theorem 2.1 implies that in full generality , the feasible set of an R-LD-LP is unlik ely to ha v e a tractable description. The situation for particular instances of the problem ma y be nicer ho wev er. Let us work through a concrete example to build some geometric in tuition. Example 2.2. Consider an R-LD-LP deﬁne d by the fol lowing data: (2.5) A =     − 1 0 0 − 1 0 1 1 1     , b =     1 1 1 3     , c =  − 1 0  , G =  0 . 6 − 0 . 4 0 . 8 0 . 5  . 4 F ollowing the con ven tion of [ 18 ], we allow paths to revisit the same no de sev eral times. 8 A. A. AHMADI AND O. G ¨ UNL ¨ UK F or k = 0 , 1 , . . . , let P k : = { x ∈ R 2 | AG k x ≤ b } , so that the fe asible solutions to the pr oblem b elong to the set S = T ∞ k =0 P k . In Figur e 2 , we show P 0 , P 1 , and P 0 ∩ P 1 . P 1 x 1 x 2 P 0 x 1 x 2 P 1 x 1 x 2 P 0 ∩ P 1 Figure 2: P 0 ∩ P 1 . In Figur e 3 , we show P 0 ∩ P 1 , P 2 , and P 0 ∩ P 1 ∩ P 2 , which happ ens to b e e qual to S . This is b e c ause ( P 0 ∩ P 1 ∩ P 2 ) ⊆ P 3 , which implies (se e L emma 2.3 in Se ction 2.1 ) that ( P 0 ∩ P 1 ∩ P 2 ) ⊆ P k for al l k > 2 . The optimal value for this example is achieve d at the rightmost vertex of S and is e qual to 1 . 1492 . x 1 x 2 P 0 ∩ P 1 x 1 x 2 P 2 P 0 ∩ P 1 ∩ P 2 x 1 x 2 Figure 3: S = P 0 ∩ P 1 ∩ P 2 . 2.1. Outer appro ximations to S . F or an integer r ≥ 0 , let (2.6) S r : = r \ k =0 { x ∈ R n | AG k x ≤ b } . In view of the deﬁnition of S in ( 2.3 ), we hav e S ⊆ . . . S r +1 ⊆ S r ⊆ . . . ⊆ S 2 ⊆ S 1 ⊆ S 0 = P . ROBUST-TO-D YNAMICS OPTIMIZA TION 9 In other w ords, the sets S r pro vide p olyhedral outer appro ximations to the feasible set S of an R-LD-LP , which get tigh ter as r increases. It can b e shown that the sets S r con v erge 5 to S in the limit. By solving linear programs that minimize c T x o v er S r , one obtains a nondecreasing sequence of low er b ounds on the optimal v alue of our R-LD-LP . Moreov er, one can show that if S r is b ounded for some r , this sequence con v erges to the optimal v alue of the R-LD-LP as r → ∞ . 6 W e hav e already seen an example where conv ergence of S r to S is ﬁnite (Example 2.2 ), and one where it is not (cf. Theorem 2.1 , part (ii)). Our goal in this subsection is to study the following t w o questions: (i) What are the barriers that can preven t con v ergence of S r to S from b eing ﬁnite (i.e., mak e S r \ S  = ∅ for all p ositiv e integers r )? (ii) When con v ergence of S r to S is ﬁnite, i.e., when S = S r ∗ for some p ositive integer r ∗ , can we pro vide an eﬃciently computable upp er b ound on r ∗ ? In regards to question (i), we show that there are three separate barriers to ﬁnite conv er- gence (meaning that in the presence of any of these conditions conv ergence ma y or may not b e ﬁnite): the matrix G having sp ectral radius larger or equal to 1 (Prop osition 2.4 and The- orem 2.1 , part (ii)), the origin b eing on the b oundary of the p olyhedron P (Prop osition 2.5 ), and the p olyhedron P b eing unbounded (Prop osition 2.6 ). In regards to question (ii), w e sho w that once these three barriers are remov ed, then S r reac hes S in a n um b er of steps that is not only ﬁnite, but upp er b ounded b y a quantit y that can b e computed in p olynomial time (Theorem 2.7 ). Before we prov e these results, let us start with a simple lemma that allows us to detect ﬁnite termination. Lemma 2.3. If S r = S r +1 for some inte ger r ≥ 0 , then S = S r . F urthermor e, for any ﬁxe d r ≥ 0 , the c ondition S r = S r +1 c an b e che cke d in p olynomial time. Pr o of. W e ﬁrst observe that condition S r = S r +1 implies that the set S r is in v ariant. If not, there w ould exist an x ∈ S r with Gx / ∈ S r . But this implies that x / ∈ S r +1 , which is a con tradiction. Inv ariance of S r implies that S r = S and the ﬁrst part of the claim is prov en. T o show the second part of the claim, note that S r +1 ⊆ S r b y deﬁnition, so only the rev erse inclusion needs to b e chec k ed. F or i ∈ { 1 , . . . , m } , let a i denote the transp ose of the i -th row of A . W e can then solv e, in p olynomial time, m linear programs (2.7) max x ∈ R n { a T i G r +1 x : x ∈ S r } , and declare that S r = S r +1 if and only if the optimal v alue of the i -th program is less than or equal to b i for all i ∈ { 1 , . . . , m } . 5 By con v ergence of the sequence { S r } to S , w e mean (see [ 46 ]) that ∀ x ∈ R n , lim r →∞ dist( x, S r ) = dist( x, S ) . Here, for a p oint x ∈ R n and a set Ω ⊆ R n , dist( x, Ω) is deﬁned as inf y ∈ Ω || y − x || . 6 The assumption on b oundedness of S r for some r cannot b e dropp ed as one can verify b y considering the R-LD-LP instance given b y A = [1 − 1] , b = 1 , c = [ − 1 0] T , G =  1 0 0 1 / 2  . Note also that when P is b ounded, this assumption is clearly satisﬁed. 10 A. A. AHMADI AND O. G ¨ UNL ¨ UK In view of this lemma, the reader can see that in Example 2.2 , the observ ation that S 2 equals S 3 allo w ed us to conclude that S equals S 2 . W e no w characterize scenarios where con v ergence of S r to S is not ﬁnite. Note that this is equiv alen t to having S r +1 ⊂ S r for all r ≥ 0 . Recall that the sp e ctr al r adius ρ ( G ) of an n × n matrix G is given b y ρ ( G ) = max {| λ 1 | , . . . , | λ n |} , where λ 1 , . . . , λ n are the (real or complex) eigen v alues of G . Theorem 2.1 , part (ii) has already sho wn that when ρ ( G ) = 1 , conv ergence of S r to S ma y not b e ﬁnite. The following simple construction sho ws that the same phenomenon can o ccur when ρ ( G ) > 1 , ev en when the set S is p olyhedral. S 0 S 1 S 4 S Figure 4: The construction in the pro of of Prop osition 2.4 . Prop osition 2.4. If ρ ( G ) > 1 , then c onver genc e of S r to S may not b e ﬁnite even when P is a b ounde d p olyhe dr on that c ontains the origin in its interior. Pr o of. Let a > 0 and consider an instance of an R-LD-LP with P = [ − 1 , 1] 2 and G =  a 0 0 1 /a  . It is easy to see that for r ≥ 0, S r = { x ∈ R 2 | − a − r ≤ x 1 ≤ a − r , − 1 ≤ x 2 ≤ 1 } . Hence, the set S is the line segment joining the p oints (0 , − 1) and (0 , 1), and conv ergence of S r to S is not ﬁnite. See Figure 4 for an illustration. Prop osition 2.5. If the origin is not c ontaine d in the interior of P , then c onver genc e of S r to S may not b e ﬁnite even when ρ ( G ) < 1 and P is b ounde d. Pr o of. Consider an instance of an R-LD-LP in R 2 with P = [0 , 1] 2 and (2.8) G = 1 2  2 / 3 − 1 / 3 − 1 / 3 2 / 3  . Note that ρ ( G ) = 1 / 2 and the origin is contained in P , but not in the interior of P . It can b e c hec k ed that (2.9) G k = 1 2 k +1 " 1 − 1 − 1 1 # + 1 3 k " 1 1 1 1 #! ROBUST-TO-D YNAMICS OPTIMIZA TION 11 for any in teger k ≥ 1 . It follo ws that G k x ≥ 0 ⇐ ⇒ | x 1 − x 2 | ≤ (1 / 3 k )( x 1 + x 2 ) , (2.10) and therefore { x ∈ R 2 | G k +1 x ≥ 0 } ⊂ { x ∈ R 2 | G k x ≥ 0 } for any in teger k ≥ 1. Similarly , 1 ≥ G k x ⇐ ⇒ 2 k +1 ≥ | x 1 − x 2 | + (1 / 3 k )( x 1 + x 2 ) , ∀ k ≥ 1 . (2.11) Observ e that for an y k ≥ 1 , if x ∈ P = [0 , 1] 2 , then 2 k +1 ≥ | x 1 − x 2 | + (1 / 3 k )( x 1 + x 2 ). Com bining this observ ation with ( 2.10 ) and ( 2.11 ), we get that (2.12) S r = r \ k =0 { x ∈ R n | G k x ∈ P } = P ∩ { x ∈ R n | | x 1 − x 2 | ≤ (1 / 3 r )( x 1 + x 2 ) } , and therefore S r strictly contains S r +1 for all r ≥ 1. This shows that con v ergence of S r to S = { x ∈ R 2 | 0 ≤ x 1 ≤ 1 , x 2 = x 1 } cannot b e ﬁnite. Figure 5 demonstrates this asymptotic con v ergence. S 0 x 1 x 2 1 1 1 S 1 x 1 x 2 1 1 S 2 x 1 x 2 1 1 S 3 x 1 x 2 1 1 Figure 5: The construction in the pro of of Prop osition 2.5 . Prop osition 2.6. If P is unb ounde d, then c onver genc e of S r to S may not b e ﬁnite even if ρ ( G ) < 1 and the origin is in the interior of P . Pr o of. Consider an instance of R-LD-LP in R 2 with P = { ( x 1 , x 2 ) ∈ R 2 | x 1 ≥ − 1 } and G as in ( 2.8 ). Note that P is a half-space that contains the origin in its interior and ρ ( G ) = 1 / 2 . Recall our notation P k = { x ∈ R n | G k x ∈ P } and S r = ∩ r k =0 P k . Our goal is to pro vide a sequence of p oin ts { z r } r ≥ 1 suc h that z r ∈ P k , ∀ k = 0 , . . . , r, but z r / ∈ P r +1 . This w ould imply that S r +1 is a strict subset of S r for all r ≥ 1, which sho ws that conv ergence of S r to S cannot b e ﬁnite. W e ﬁrst start by c haracterizing the sets P k for k ≥ 1 . In view of ( 2.9 ), we ha v e G k  x 1 x 2  = 1 2 k +1  x 1 (1 + 1 3 k ) + x 2 ( 1 3 k − 1) x 1 ( − 1 + 1 3 k ) + x 2 ( 1 3 k + 1)  . Hence, it is easy to chec k that for k ≥ 1, P k is the follo wing half-space: P k = { ( x 1 , x 2 ) ∈ R 2 | x 2 ≤ l k x 1 + i k } , where l k := 3 k + 1 3 k − 1 and i k := 2 · 6 k 3 k − 1 . 12 A. A. AHMADI AND O. G ¨ UNL ¨ UK W e now deﬁne, for r ≥ 1, the co ordinates of the p oints z r = ( z r, 1 , z r, 2 ) T to b e z r, 1 = i r +2 − i r l r − l r +2 and z r, 2 = i r +2 + l r +2 ·  i r +2 − i r l r − l r +2  . Fix r ≥ 1. Note that z r, 1 = 2 r − 3 (3 r +3 − 35) ≥ 1 4 (3 4 − 35) ≥ − 1 , which implies that z r ∈ P . T o sho w that z r ∈ P k , ∀ k = 1 , . . . , r , w e need to show that z r, 2 ≤ l k z r, 1 + i k , for k = 1 , . . . , r. This is the same as showing that i r +2 − i r l r − l r +2 ≥ i r +2 − i k l k − l r +2 , for k = 1 , . . . , r, whic h is equiv alen t to sho wing that the diﬀerence of the tw o ratios, i.e., (3 r +2 − 1)(2 k +3 3 k + 2 r 3 r +3 − 2 r +3 3 k − 2 r 3 k +3 ) 2 3 (3 r +2 − 3 k ) (2.13) is nonnegative for k = 1 , . . . , r . As we hav e 2 k +3 3 k + 2 r 3 r +3 − 2 r +3 3 k − 2 r 3 k +3 = 2 r 3 k (2 3 − ( r − k ) + 3 3+( r − k ) − 2 3 − 3 3 ) and as x ≥ 0 ⇒ 2 3 − x + 3 3+ x − 2 3 − 3 3 ≥ 0, w e hav e z r ∈ P k , ∀ k = 1 , . . . , r. T o show that z r / ∈ P r +1 , we simply need to sho w that i r +2 − i r l r − l r +2 < i r +2 − i r +1 l r +1 − l r +2 . Replacing k b y r + 1 in ( 2.13 ), this is equiv alen t to sho wing that 2 r +4 3 r +1 + 2 r 3 r +3 − 2 r +3 3 r +1 − 2 r 3 r +4 < 0 . As 2 r +4 3 r +1 + 2 r 3 r +3 − 2 r +3 3 r +1 − 2 r 3 r +4 = 2 r 3 r +1 (2 4 + 3 2 − 2 3 − 3 3 ) = − 10 · 2 r 3 r +1 , this inequality clearly holds. S 0 x 1 x 2 S 1 x 1 x 2 S 2 x 1 x 2 S 3 x 1 x 2 Figure 6: The construction in the pro of of Prop osition 2.6 . ROBUST-TO-D YNAMICS OPTIMIZA TION 13 2.1.1. A p olynomially-computable upper b ound on the numb er of steps to conver- gence. In this subsection, we show that when (i) ρ ( G ) < 1, (ii) P is b ounded, and (iii) the origin is in the interior of P , then S = S r for some integer r that can b e computed in time p olynomial in the size of the R-LD-LP data. Note that in view of Prop ositions 2.4 , 2.5 , 2.6 , when one of the assumptions (i), (ii), (iii) ab o v e do es not hold, one can ﬁnd examples where con v ergence of S r to S is not ev en ﬁnite. Arguably , the condition ρ ( G ) < 1 accounts for one of the more interesting settings of an R-LD-LP . Indeed, if ρ ( G ) > 1 , then tra jectories of the dynamical system x k +1 = Gx k starting from all but a measure zero set of initial conditions go to inﬁnity and hence, at least when the p olyhedron P is b ounded, the feasible set S of our R-LD-LP can nev er b e full dimensional. The b oundary case ρ ( G ) = 1 is more delicate. Here, the system tra jectories can stay b ounded or go to inﬁnit y depending on the geometric/algebraic multiplicit y of the eigenv alues of G with absolute v alue equal to one. Even in the b ounded case, we hav e sho wn already in Theorem 2.1 , part (ii) that the feasible set of an R-LD-LP may not b e p olyhedral. Hence the optimal v alue of an R-LD-LP may not ev en b e a rational num b er (consider, e.g., the set S asso ciated with Figure 1 with c = (1 , 1) T ). Note also that when ρ ( G ) < 1, we hav e lim k →∞ G k x 0 = 0 for all x 0 ∈ R n and hence if the origin is not in P , then the feasible set of the R-LD-LP is empty . As a consequence, the assumption that the origin b e in P is reasonable and that it b e in the in terior of P is only sligh tly stronger (and cannot b e a v oided). F or the conv enience of the reader, we next give the standard deﬁnitions (see, e.g., [ 47 ]) on sizes of rational data that we will use in the follo wing result. W e say that the size of a rational n um b er r = p/q where p, q ∈ Z (and are relativ ely prime), is 1 + ⌈ log 2 ( | p | + 1) ⌉ + ⌈ log 2 ( | q | + 1) ⌉ . W e denote the size of r by σ ( r ) and note that 1 / | r | , | r | ≤ 2 σ ( r ) . Similarly , the size of a rational vector (or matrix) is deﬁned to b e the sum of the sizes of its comp onents plus the pro duct of its dimensions. It is w ell known that multiplying tw o matrices gives a matrix of size p olynomially b ounded b y the sizes of the initial matrices. The in v erse of a nonsingular rational matrix has size polynomially b ounded b y the size of the matrix. Similarly , an y system of rational linear equations has a solution of size p olynomially b ounded b y the size of the data deﬁning the equality system [ 47 ]. In addition, this solution can b e computed in p olynomial time. Consequen tly , if a linear program deﬁned by rational data has an optimal solution, then it has one with size p olynomial in the data deﬁning the LP . Clearly the optimal v alue of the LP has size p olynomial in the data as well. Finally , w e remark that giv en a p olyhedron P = { x ∈ R n | Ax ≤ b } , one can chec k whether P contains the origin in its interior and whether P is b ounded in time p olynomial in σ ( A, b ). The former task simply requires c hec king if the en tries of b are all positive, and the latter can b e carried out e.g. by minimizing and maximiz ing each co ordinate x i o v er P and chec king if the optimal v alues of the resulting LPs are all ﬁnite. One can also chec k if the sp ectral radius of a square matrix G is less than one in time p olynomial in σ ( G ) [ 19 , Section 2.6]. Theo rem 2.7. L et σ ( A, b, G ) denote the size of A , b, and G . L et S and S r b e as in ( 2.3 ) and ( 2.6 ) r esp e ctively. If ρ ( G ) < 1 , P = { x ∈ R n | Ax ≤ b } is b ounde d, and the origin is in the interior of P , then S = S r for some nonne gative inte ger r that c an b e c ompute d in time p olynomial in σ ( A, b, G ) . F urthermor e, for any ﬁxe d r ational numb er ρ ∗ < 1 , R-LD-LPs with 14 A. A. AHMADI AND O. G ¨ UNL ¨ UK ρ ( G ) ≤ ρ ∗ c an b e solve d in time p olynomial in σ ( A, b, c, G, ρ ∗ ) . 7 Pr o of. F or the ﬁrst claim, w e prov e that r can b e computed in p olynomial time by paying atten tion to the bit-size complexit y of the n um b ers in volv ed in each of the follo wing four steps: • First, w e ﬁnd a p ositive deﬁnite matrix M ≻ 0 8 and ellipsoids E ( α ) = { x ∈ R n | x T M x ≤ α } that are in v arian t under the dynamics x k +1 = Gx k . • Next, we ﬁnd scalars α 1 , α 2 > 0 suc h that E ( α 1 ) ⊆ P ⊆ E ( α 2 ). • Then we compute a “shrink age factor” γ ∈ (0 , 1), whic h gives a low er b ound on the amoun t our ellipsoids E ( α ) shrink under one application of G . • Finally , using α 1 , α 2 , γ , w e compute a nonnegativ e integer r such that G r E ( α 2 ) ⊆ E ( α 1 ) . This will imply that S = S r . The second claim of the theorem then follows by adapting Steps 1 and 3 exploiting the fact that ρ ∗ < 1 is ﬁxed. Step 1. Computing an inv ariant ellipsoid. T o ﬁnd an inv ariant ellipsoid for the dynam- ical system x k +1 = Gx k , we solv e the well-kno wn Lyapuno v equation (2.14) G T M G − M = − I , for the symmetric matrix M , where I here is the n × n identit y matrix. Since ρ ( G ) < 1, this linear system is guaranteed to hav e a unique solution (see e.g. [ 22 , Chap. 14]). Note that M is rational, can b e computed in p olynomial time, and has size p olynomially b ounded b y σ ( G ) . F urther, we claim that M m ust b e p ositive deﬁnite. T o see this, supp ose w e had y T M y ≤ 0, for some y ∈ R n , y  = 0 . Multiplying ( 2.14 ) from left and righ t b y y T and y , w e see that y T G T M Gy ≤ − y T y < 0. In fact, y T ( G k ) T M G k y ≤ − y T y < 0, for all k ≥ 1. But since ρ ( G ) < 1 , we must ha v e G k y → 0 as k → ∞ , and hence y T ( G k ) T M G k y → 0 as k → ∞ , a contradiction. Since M ≻ 0, the sets E ( α ) : = { x ∈ R n | x T M x ≤ α } are b ounded for all α ≥ 0 . F urthermore, b ecause G T M G ≺ M in view of ( 2.14 ), if x ∈ E ( α ) for some α ≥ 0, then Gx ∈ E ( α ) , and hence G k x ∈ E ( α ) for all in tegers k ≥ 0. Step 2. Computing inner and outer ellipsoids. W e next compute scalars α 2 > α 1 > 0 suc h that the ellipsoids E ( α 1 ) = { x ∈ R n | x T M x ≤ α 1 } and E ( α 2 ) = { x ∈ R n | x T M x ≤ α 2 } satisfy E ( α 1 ) ⊆ P ⊆ E ( α 2 ) . F or i = 1 , . . . , m , let a T i x ≤ b i denote the i -th deﬁning inequality of P . F or eac h i , compute a scalar η i as the optimal v alue of the following conv ex program: η i : = max x ∈ R n { a T i x : x T M x ≤ 1 } . 7 Our preliminary version of this work [ 2 ] unfortunately did not hav e the assumption that ρ ∗ b e ﬁxed, which is needed in our pro of. The statement in the pro of of Theorem 3.1 of [ 2 ] that the integer r computed there has p olynomial size is correct, but that linear optimization ov er S r is an LP of p olynomial size is incorrect. W e ask that the reader refer to the statemen t and pro of of Theorem 2.7 of the current pap er instead. 8 W e use the standard notation A ≻ 0 to denote that a matrix A is p ositiv e deﬁnite (i.e., has positive eigen v alues), and A ⪰ 0 to denote that A is p ositiv e semideﬁnite (i.e., has nonnegative eigenv alues). ROBUST-TO-D YNAMICS OPTIMIZA TION 15 The optimal v alue here can b e computed in closed form: η i = q a T i M − 1 a i . Note that M − 1 exists as M ≻ 0. W e then let (2.15) α 1 = min i ∈{ 1 ,...,m }  b 2 i /η 2 i  = min i ∈{ 1 ,...,m }  b 2 i a T i M − 1 a i  . This ensures that ∀ i ∈ { 1 , . . . , m } , ∀ x ∈ E ( α 1 ), w e hav e a T i x ≤ b i . Hence, E ( α 1 ) ⊆ P . Note that α 1 > 0 since M − 1 ≻ 0 , and b i > 0 for i = 1 , . . . , m as the origin is in the in terior of P . Moreov er, the size of α 1 is p olynomially bounded by σ ( A, b, G ) as the size of the num b ers used in the calculation are all p olynomially b ounded b y σ ( A, b, G ). W e next compute a scalar α 2 > 0 suc h that P ⊆ E ( α 2 ) . As P is a p olytope, P ⊆ { x ∈ R n | l i ≤ x i ≤ u i } , where the scalars u i , l i ∈ R , for i = 1 , . . . , n , can b e obtained by solving LPs that minimize/ maximize eac h co ordinate x i o v er P . Note that this can b e done in polynomial time. W e can no w b ound x T M x = P i,j M i,j x i x j term by term to obtain (2.16) α 2 = n X i =1 n X j =1 max {| M i,j u i u j | , | M i,j l i l j | , | M i,j u i l j | , | M i,j l i u j |} . Clearly for any x ∈ P , we hav e x T M x ≤ α 2 and therefore P ⊆ E ( α 2 ). As optimal v alues of linear programs are b ounded p olynomially by the size of their data, we hav e the size of α 2 p olynomially b ounded by σ ( A, b, G ). Step 3. Computing a shrink age factor. Next w e argue that the num b er (2.17) γ = 1 − 1 max i ∈{ 1 ,...,n } { M ii + P j  = i | M i,j |} satisﬁes G T M G ⪯ γ M . Observ e that for any x ∈ R n , the Ly apuno v equation ( 2.14 ) implies x T G T M Gx = x T M x − x T x ≤ (1 − η ) x T M x, pro vided that η > 0 is a scalar that satisﬁes the inequalit y η x T M x ≤ x T x for all x ∈ R n . Let λ max ( M ) denote the largest eigenv alue of the matrix M . Note that any η ≤ 1 /λ max ( M ) satisﬁes the ab o v e inequality . By Gershgorin’s circle theorem, w e hav e λ max ( M ) ≤ max i ∈{ 1 ,...,n } { M ii + X j  = i | M i,j |} . Using this upp er b ound, we establish that G T M G ⪯ γ M for γ as given in ( 2.17 ). 16 A. A. AHMADI AND O. G ¨ UNL ¨ UK W e observe that as M ≻ 0, and hence G T M G ⪰ 0 , w e hav e M − I = G T M G ⪰ 0. This implies that M ii ≥ 1 for all i ∈ { 1 , . . . , n } and therefore γ is indeed a n um b er in [0 , 1). Note also that the size of γ is p olynomially b ounded by σ ( G ) . Step 4. Computing the n um b er of steps to con v ergence. Using α 1 , α 2 and γ , we no w compute an in teger ¯ r such that γ ¯ r E ( α 2 ) ⊆ E ( α 1 ) and therefore all points inside the outer ellipsoid E ( α 2 ) are guaran teed to b e within the inner ellipsoid E ( α 1 ) after at most ¯ r steps. As G T M G ⪯ γ M , w e hav e ( G k ) T M G k ⪯ γ k M . Hence, if x ∈ E ( α 2 ), i.e., x T M x ≤ α 2 , then x T ( G k ) T M G k x ≤ γ k α 2 . Clearly γ ¯ r α 2 ≤ α 1 for ¯ r =  log( α 1 /α 2 ) log γ  =  log( α 2 /α 1 ) log(1 /γ )  ≤  ( α 2 /α 1 ) − 1 (1 − γ )  = : r , where the inequalit y ab o v e uses the fact that 1 − (1 /a ) ≤ log a ≤ a − 1 for an y scalar a > 0. Therefore, an y p oint x ∈ E ( α 2 ) satisﬁes G r x ∈ E ( α 1 ). As E ( α 1 ) is inv ariant and E ( α 1 ) ⊆ P , w e conclude that if G t x ∈ P for t = 1 , . . . , r , then G t x ∈ P for all t ≥ r . This establishes that S = S r (in fact, we hav e sho wn that S = S r = S ¯ r ). Note that the n um b ers α 1 , α 2 , γ , and hence r, can b e computed in time p olynomial in σ ( A, b, G ) . This completes the ﬁrst part of the pro of. Solving R-LD-LP in p olynomial time. W e next sho w that the num ber of steps to con v ergence is itself p olynomially b ounded by σ ( A, b, G, ρ ∗ ) when ρ ( G ) is upp er b ounded by a rational constant ρ ∗ < 1. This would imply that after a polynomial num b er of steps r , w e would hav e S = S r and therefore the inequalities describing S could b e written down in p olynomial time. As an y linear function c T x can b e optimized o ver a p olyhedron in p olynomial time, this w ould prov e the second claim of the theorem. T o this end, w e compute the inv arian t ellipsoid E in Step 1 and the shrink age factor γ in Step 3 slightly diﬀeren tly . T o ﬁnd an inv ariant ellipsoid for G , w e now solve the Ly apuno v equation ˆ G T ˆ M ˆ G − ˆ M = − I , where ˆ G = (1 / ˆ ρ ) G and ˆ ρ = (1 + ρ ∗ ) / 2. As ρ ( G ) < ˆ ρ < 1, w e hav e ρ ( ˆ G ) < 1 and therefore the ab o v e equation has a unique solution. Moreov er, the size of ˆ M is p olynomially b ounded b y σ ( A, b, ˆ G ) and therefore by σ ( A, b, G, ρ ∗ ). (Recall also that ρ ∗ is a constant throughout.) As ˆ G T ˆ M ˆ G − ˆ M is negativ e deﬁnite, ˆ M ≻ ˆ G T ˆ M ˆ G = 1 ˆ ρ 2 G T ˆ M G whic h readily gives the shrink age factor ˆ γ = ˆ ρ 2 . W e can now compute ˆ α 1 and ˆ α 2 using Equations ( 2.15 ) and ( 2.16 ) with ˆ M . Clearly the sizes of b oth ˆ α 1 and ˆ α 2 are p olynomially b ounded b σ ( A, b, ˆ G ) and therefore by σ ( A, b, G, ρ ∗ ). W e observe with the same argument as b efore that with r =  log( ˆ α 2 / ˆ α 1 ) log(1 / ˆ γ )  , w e must hav e S = S r . Note that the size of ˆ α 1 / ˆ α 2 , i.e., log ( ˆ α 1 / ˆ α 2 ), is p olynomially b ounded b y σ ( A, b, G, ρ ∗ ). As ˆ γ is a constant, r is indeed p olynomially b ounded by σ ( A, b, G, ρ ∗ ). ROBUST-TO-D YNAMICS OPTIMIZA TION 17 T o summarize, under the assumptions of Theorem 2.7 , we ha v e provided a pseudo- p olynomial time algorithm 9 for R-LD-LP , and a p olynomial algorithm for all instances where the sp ectral radius of G is upp er b ounded b y a constant less than one. W e end this subsection b y remarking that since a conv ex quadratic function can b e minimized o v er a polyhedron in p olynomial time [ 34 ], the same complexity guaran tees carry ov er to a generalization of R-LD-LP where the linear ob jective function is replaced b y a con v ex quadratic one. 2.2. Inner app roximations to S . In this subsection, w e fo cus on the computation of a sequence of inner appro ximations to the feasible set S that pro duce fe asible solutions to the R-LD-LP in each step. By minimizing c T x ov er these sets, one obtains upp er b ounds on the optimal v alue of the R-LD-LP . (Note that p oints b elonging to the outer approximations S r ma y not b e feasible to the R-LD-LP , unless we wait long enough for S r to coincide with S .) Motiv ated b y the analysis in Section 2.1 , we are in terested in the remainder of this section in the setting where ρ ( G ) < 1 and P is a b ounded p olyhedron that contains the origin in its interior. Some of the statemen ts in our lemmas b elow ho w ev er do not need all three assumptions. Recall the notation S r for the outer-approximating p olyhedra deﬁned in ( 2.6 ). T o ﬁnd a family of inner approximations to S , we assume that an inv ariant set E ⊆ P (with resp ect to G ) is giv en. W e will discuss the eﬃcien t computation of the set E later (cf. Section 2.2.1 and Section 2.2.2 ). Deﬁne S − 1 := R n and for an y integer r ≥ 0, let (2.18) I r ( E ) = S r − 1 ∩ { x ∈ R n | G r x ∈ E } . Note that I 0 ( E ) = E by deﬁnition. W e next argue that the sets I r ( E ) are nested and con tained in S . Lemma 2.8. L et E ⊆ P b e invariant with r esp e ct to G . Then, for any inte ger r ≥ 0 , (i) I r ( E ) ⊆ S , (ii) I r ( E ) ⊆ I r +1 ( E ) , and (iii) if I r ( E ) = I r +1 ( E ) , then I k ( E ) = I r ( E ) for al l k ≥ r . Pr o of. First note that as I 0 ( E ) = E ⊆ P , and E is inv arian t, w e ha v e I 0 ( E ) ⊆ S . When r ≥ 1, if x ∈ I r ( E ) then x ∈ S r − 1 and therefore G t x ∈ P for t = 0 , . . . , r − 1 . In addition, as the set E ⊆ P is in v ariant, G r x ∈ E implies that G r + t x ∈ E ⊆ P for t ≥ 0 . Consequently , G t x ∈ P for all t ≥ 0 and therefore I r ( E ) ⊆ S . T o see that I r ( E ) ⊆ I r +1 ( E ), note that if x ∈ I r ( E ) then x ∈ S r − 1 and G r x ∈ E ⊆ S . As G r x ∈ S implies AG r x ≤ b , we hav e x ∈ S r . F urthermore, as E is in v ariant, if G r x ∈ E then G r +1 x ∈ E as w ell. Consequently , x ∈ I r +1 ( E ) as desired. T o prov e that the last claim also holds, we will argue that if I r ( E ) = I r +1 ( E ), then I r +1 ( E ) = I r +2 ( E ). As I r +1 ( E ) ⊆ I r +2 ( E ), w e need to sho w that I r +1 ( E ) ⊇ I r +2 ( E ). Assume I r +1 ( E ) ⊇ I r +2 ( E ), and let x ∈ I r +2 ( E ) \ I r +1 ( E ). In this case, as x ∈ I r +2 ( E ), we ha v e x ∈ S r +1 and G r +2 x ∈ E . F urthermore, as x ∈ I r +1 ( E ), and x ∈ S r ⊇ S r +1 , w e also ha v e G r +1 x ∈ E . Now consider the p oin t y = Gx . Clearly , y ∈ S r as x ∈ S r +1 . F urthermore, w e ha v e G r +1 y ∈ E as G r +2 x ∈ E and therefore, y ∈ I r +1 ( E ) . Ho wev er, G r y ∈ E as 9 W e recall that a pseudo-p olynomial time algorithm is an algorithm whose running time is p olynomial in the numeric v alue of the input (the largest in teger presen t in the input), but not necessarily in the size of the input (the num b er of bits required to represent the input). 18 A. A. AHMADI AND O. G ¨ UNL ¨ UK G r +1 x ∈ E and therefore y ∈ I r ( E ) . This contradicts the assumption that I r ( E ) = I r +1 ( E ) as y ∈ I r +1 ( E ) \ I r ( E ). W e conclude that for all r ≥ 0 (2.19) I r ( E ) ⊆ I r +1 ( E ) ⊆ S ⊆ S r +1 ⊆ S r pro vided that the set E is inv ariant with resp ect to G and E ⊆ P . Also note that we did not mak e any assumptions on P or on G for the inclusion relationships in ( 2.19 ) to hold. Lemma 2.9. L et P b e a b ounde d p olyhe dr on, ρ ( G ) < 1 , and E ⊆ P b e invariant with r esp e ct to G . F urthermor e, assume that E c ontains the origin in its interior. Under these assumptions, if I r ( E ) = I r +1 ( E ) , then S = I r ( E ) . Pr o of. Let L k : = { x ∈ R n | G k x ∈ E } and recall that I k ( E ) = S k − 1 ∩ L k and therefore I k ( E ) ⊆ S k − 1 for all k > 0. Also note that lim k →∞ S k = S and as I r ( E ) = I r +1 ( E ) by assumption, Lemma 2.8 implies that lim k →∞ I k ( E ) = I r ( E ). Therefore, taking the limit, w e obtain: I r ( E ) = S ∩ lim k →∞ L k . Consequen tly , to prov e the claim, we need to argue that lim k →∞ L k ⊇ S . W e will actually sho w that lim k →∞ L k ⊇ P , whic h is suﬃcient as P ⊇ S . As ρ ( G ) < 1 , follo wing the steps in the pro of of Theorem 2.7 , w e can ﬁnd a positive deﬁnite matrix M suc h that the ellipsoid E ( β ) = { x ∈ R n | x T M x ≤ β } is inv ariant under the linear dynamics G for all β > 0. As E is full-dimensional and con tains the origin is its in terior, there exists a scalar α 1 > 0 such that E ⊇ E ( α 1 ). In addition, w e can compute a scalar α 2 > 0 suc h that E ( α 2 ) ⊇ P . Therefore, w e hav e E ( α 1 ) ⊆ E ⊆ P ⊆ E ( α 2 ) . F urthermore, the shrink age factor giv en by equation ( 2.17 ) implies that for some nonnegativ e in teger m, all x ∈ E ( α 2 ) satisfy G k x ∈ E ( α 1 ) for all k ≥ m . As E ( α 1 ) ⊆ E and P ⊆ E ( α 2 ) , this implies that if x ∈ P , then G m x ∈ E and therefore L k ⊇ P for all k ≥ m . Consequently , lim k →∞ L k ⊇ P , and S = I r ( E ) as desired. The pro of of Lemma 2.9 sho ws that for any in v ariant set E ⊆ P, one can compute a nonnegativ e in teger m E suc h that G m E x ∈ E for all x ∈ P . This implies that L k ⊇ S for all k ≥ m E . In addition, Theorem 2.7 shows that for some nonnegativ e in teger r , w e ha v e S = S k , ∀ k ≥ r . Consequently , for all k ≥ max { m E , r } we ha v e L k ⊇ S = S k , which we formally state next. Co rolla ry 2.10. Under the assumptions of L emma 2.9 , ther e is a nonne gative inte ger t such that S = I r ( E ) for al l r ≥ t . 2.2.1. Computation of I r ( E ) . The construction of the sets I r ( E ) requires access to an in v ariant set E ⊆ P . F or an R-LD-LP with ρ ( G ) < 1, an inv ariant set for the dynamics that is alwa ys guaranteed to exist is an ellipsoid E = { x ∈ R n | x T M x ≤ α } . T o ﬁnd the p ositive deﬁnite matrix M (that ensures G T M G ⪯ M ) and the p ositiv e scalar α (that ensures E ⊆ P ), ROBUST-TO-D YNAMICS OPTIMIZA TION 19 one can follow the metho dology describ ed in steps 1 and 2 of the pro of of Theorem 2.7 . Note that these t w o steps only inv olve matrix inv ersion and basic arithmetic op erations. With M and α ﬁxed, one can solve the following sequence of con v ex quadratic programs, minimize x ∈ R n c T x s.t. ( G r x ) T M ( G r x ) ≤ α, (2.20) AG k x ≤ 1 , k = 0 , . . . , r − 1 , indexed by an integer 10 r ≥ 0. The feasible sets of these optimization problems are the sets I r ( E ) as deﬁned in ( 2.18 ), which reside inside the feasible set S of our R-LD-LP . Hence, the optimal v alues of these conv ex quadratic programs are upp er b ounds on the optimal v alue of the R-LD-LP . By Lemma 2.8 , these upp er b ounds monotonically improv e with r , and by Corollary 2.10 , they reac h the optimal v alue of the R-LD-LP in a ﬁnite num b er of steps. Although this approach is simple and conv ergent, it is sub optimal in terms of the quality of the upp er b ounds that it returns in each iteration. W e explain how one can do b etter next. 2.2.2. Computation of imp roved inner app ro ximations. Our impro vemen t ov er the algo- rithm suggested in the last subsection is based on answ ers to the following tw o basic questions: (i) Instead of ﬁnding any inv ariant ellipsoid E and then optimizing o v er the sets I r ( E ) gener- ated by it, can we search for an “optimal” in v arian t ellipsoid at the same time as we optimize o v er I r ( E )? (ii) Instead of w orking with a ﬁxed inv arian t ellipsoid throughout the hierarch y , can we reoptimize the ellipsoid at each iteration? W e show here that semideﬁnite pr o gr amming (SDP) can ac hiev e b oth of these goals at once. Let r = 0 , 1 , . . . b e the index of our hierarch y . At step r , the strategy is to ﬁnd an ellipsoid E r , deﬁned as the unit sublevel set of a quadratic form x T H r x , which satisﬁes the follo wing prop erties: 1. The set E r is inv ariant under the dynamics x k +1 = Gx k . 2. The set E r is contained in the p olytop e P . 3. Among all ellipsoids that ha v e the previous tw o prop erties, E r is one that gives the minim um v alue of c T x as x ranges o v er the p oints in R n that land in E r after r steps and do not lea ve P b efore doing so. (The set of suc h points will b e denoted b y I r ( E r ).) As w e are under the running assumption that the origin is in the in terior of our p olytop e P , the vector b ∈ R m in the description { x ∈ R n | Ax ≤ b } of P is element wise p ositiv e. Hence, b y rescaling, w e can without loss of generality take b to b e the all ones vector. With this in mind, here is a mathematical description of the ab ov e optimization problem 11 : 10 Note that when r = 0, the ﬁnal set of constraints drop out. 11 W e use the notation S n × n to refer to the set of n × n real symmetric matrices. 20 A. A. AHMADI AND O. G ¨ UNL ¨ UK minimize x ∈ R n ,H ∈ S n × n c T x s.t. H ≻ 0 , (2.21) G T H G ⪯ H , ∀ z ∈ R n , z T H z ≤ 1 = ⇒ Az ≤ 1 , ( G r x ) T H ( G r x ) ≤ 1 , AG k x ≤ 1 , k = 0 , . . . , r − 1 . If the pair ( x r , H r ) is an optimal solution to this problem, then w e let E r = { z ∈ R n | z T H r z ≤ 1 } , (2.22) I r ( E r ) = { z ∈ R n | ( G r z ) T H r ( G r z ) ≤ 1 , AG k z ≤ 1 , k = 0 , . . . , r − 1 } , and x r will be our candidate sub optimal solution to R-LD-LP . There are t w o c hallenges to o v ercome with the formulation in ( 2.21 ). First, the constraint ∀ z ∈ R n , z T H z ≤ 1 = ⇒ Az ≤ 1 needs to b e rewritten to remov e the universal quan tiﬁer. Second, the decision v ariables x and H are m ultiplying each other in the constraint ( G r x ) T H ( G r x ) ≤ 1, which again mak es the constraint nonconv ex. Nevertheless, w e show next that one can get around these issues and formulate problem ( 2.21 ) exactly as an SDP . The main ingredients of the pro of are Sch ur complemen ts, p olar dualit y theory of conv ex sets (see e.g. [ 8 ], [ 44 ]), and dualit y of linear dynamical systems under transp osition of the matrix G . F urthermore, we establish that the feasible solutions to R-LD-LP that are pro duced by our SDPs b ecome optimal in a num b er of steps that can b e computed in p olynomial time. Theo rem 2.11. Supp ose ρ ( G ) < 1 and the set P = { x ∈ R n | Ax ≤ 1 } is b ounde d. L et a i denote the tr ansp ose of the i -th r ow of the matrix A ∈ R m × n and c onsider the fol lowing semideﬁnite pr o gr am: minimize x ∈ R n ,Q ∈ S n × n c T x s.t. Q ≻ 0 , (2.23) GQG T ⪯ Q, a T i Qa i ≤ 1 , i = 1 , . . . , m,  Q G r x ( G r x ) T 1  ⪰ 0 , AG k x ≤ 1 , k = 0 , . . . , r − 1 . Then, ROBUST-TO-D YNAMICS OPTIMIZA TION 21 (i) the optimal values of pr oblems ( 2.21 ) and ( 2.23 ) ar e the same, the optimal ve ctors x r in the two pr oblems ar e the same, and the optimal matric es H r and Q r ar e r elate d via Q r = H − 1 r . (ii) the optimal values of the SDPs in ( 2.23 ) pr ovide upp er b ounds on the optimal value of the R-LD-LP, ar e nonincr e asing with r , and r e ach the optimal value of the R-LD-LP in a (ﬁnite) numb er of steps ¯ r which c an b e c ompute d in time p olynomial in σ ( A, G ) . Mor e over, any optimal solution x r to the SDP in ( 2.23 ) with r ≥ ¯ r is an optimal solution to R-LD-LP. Pr o of. (i) W e show that a pair ( x, H ) is feasible to ( 2.21 ) if and only if the pair ( x, H − 1 ) is feasible to ( 2.23 ). Indeed, w e ha v e H ≻ 0 ⇔ H − 1 ≻ 0 as the eigen v alues of H − 1 are the in v erse of the eigenv alues of H . Moreo v er, b y tw o applications of the Sch ur complement (see, e.g., [ 20 , App endix A.5.5]), we observ e that G T H G ⪯ H ⇔  H − 1 G G T H  ⪰ 0 ⇔ GH − 1 G T ⪯ H − 1 . W e also ha v e that ( G r x ) T H ( G r x ) ≤ 1 ⇔  H − 1 G r x ( G r x ) T 1  ⪰ 0 , due to the Sch ur complement once again. Recall no w that for a set T ⊆ R n , its p olar dual T ◦ is deﬁned as T ◦ : = { y ∈ R n | y T x ≤ 1 , ∀ x ∈ T } . Let E : = { z ∈ R n | z T H z ≤ 1 } and P : = { z ∈ R n | Az ≤ 1 } . One can v erify that (i) E ⊆ P ⇔ P ◦ ⊆ E ◦ , (ii) E ◦ = { y ∈ R n | y T H − 1 y ≤ 1 } , and (iii) P ◦ = con v { a 1 , . . . , a m } , where conv here denotes the con v ex hull op eration. Hence we hav e ( ∀ z ∈ R n , z T H z ≤ 1 = ⇒ Az ≤ 1) ⇔ a T i H − 1 a i ≤ 1 , i = 1 , . . . , m. (ii) The statemen t that the optimal v alue of ( 2.23 ) is an upper b ound on the optimal v alue of the R-LD-LP follows from the fact that this SDP is constraining the optimal solution x r to b e in I r ( E r ), as deﬁned in ( 2.22 ), which is contained in S b y construction (cf. Lemma 2.8 ). F urthermore, if a pair ( x, Q ) is feasible to the SDP in ( 2.23 ) at level r , then it is also feasible to the SDP at level r + 1. This is b ecause E : = { y ∈ R n | y T Q − 1 y ≤ 1 } ⊆ P , and G r x ∈ E ⇒ G r +1 x ∈ E by inv ariance of E under G . Hence the claim ab out the monotonic improv ement of the upp er b ounds follows. T o pro v e the statement ab out ﬁnite termination of this SDP hierarch y in a p olynomially- computable num b er of steps, let M ≻ 0 b e the unique solution to the linear s ystem G T M G − M = − I , α 1 > 0 b e as in ( 2.15 ) with b 1 = · · · = b m = 1, α 2 b e as in ( 2.16 ), γ b e as in ( 2.17 ), and ¯ r =  ( α 2 /α 1 ) − 1 (1 − γ )  . The pro of of Theorem 2.7 already shows that this num ber can b e computed in p olynomial time. Let ¯ x b e any optimal solution to the R-LD-LP . W e claim that the pair ( ¯ x, α 1 M − 1 ) is a feasible solution to the SDP in ( 2.23 ) with r = ¯ r . Clearly , the constraints AG k ¯ x ≤ 1 , k = 0 , . . . , ¯ r − 1 are satisﬁed as ¯ x ∈ S . Moreov er, the pro of of Theorem 2.7 sho ws that the set ¯ E : = { y ∈ R n | y T M α 1 y ≤ 1 } is contained in P and is such that G ¯ r x ∈ ¯ E , ∀ x ∈ P . This, together 22 A. A. AHMADI AND O. G ¨ UNL ¨ UK with the equation G T M G − M = − I and the fact that ¯ x ∈ P , implies that the pair ( ¯ x, M α 1 ) is feasible for the problem in ( 2.21 ) with r = ¯ r . In view of the pro of of part (i) of the current theorem, our claim about feasibilit y of ( ¯ x, α 1 M − 1 ) to the SDP in ( 2.23 ) with r = ¯ r follo ws. T o ﬁnish the pro of, let ( x ¯ r , Q ¯ r ) b e an optimal solution to this SDP . W e m ust hav e c T x ¯ r ≤ c T ¯ x as w e hav e just argued ¯ x is feasible to the SDP . Y et c T x ¯ r ≥ c T ¯ x as I ¯ r ( E ¯ r ) ⊆ S and x ¯ r ∈ I ¯ r ( E ¯ r ). Hence, the optimal v alue of the SDP matches the optimal v alue of the R-LD-LP for all r ≥ ¯ r . Consequen tly , optimal solutions x r to the SDP m ust b e optimal to the R-LD-LP for all r ≥ ¯ r as they ac hiev e the optimal v alue and b elong to S . W e observ e that the size of the semideﬁnite constrain ts in ( 2.23 ), whic h are the most exp ensiv e constraints in that optimization problem, does not grow with r . Let us now give an example. Example 2.12. Consider an R-LD-LP deﬁne d by the fol lowing data: A =       1 0 − 1 . 5 0 0 1 0 − 1 1 1       , b =     1 1 1 1     , c =  − 0 . 5 − 1  , G = 4 5  cos( θ ) sin( θ ) , − sin( θ ) cos( θ )  wher e θ = π 6 . In Figur e 7 , we plot the inner appr oximations I r and outer appr oximations S r to S for r = 0 (on the left) and r = 1 (on the right). Note that when r = 0 , S 0 is simply P . We also plot the optimal solution to the pr oblem of minimizing c T x over S 0 (r esp. I 0 ) in Figur e 7a and over S 1 (r esp. I 1 ) in Figur e 7b . We r emark that the solutions do not c oincide for r = 0 but they do for r = 1 . Henc e, our metho d c onver ges in one step. This is further evidenc e d by the se quenc e of lower and upp er b ounds on the optimal value of the R-LD-LP given in T able 2 , which shows that we have r e ache d the exact optimal value at r = 1 . (a) r = 0 (b) r = 1 Figure 7: Outer and inner approximations to the feasible set of the R-LD-LP in Example 2.12 . ROBUST-TO-D YNAMICS OPTIMIZA TION 23 (a) The set I 1 ( E 1 ) (b) The set E 1 Figure 8: The set I 1 ( E 1 ) is the set of p oin ts in P that land in the ellipsoid E 1 after one application of G . r = 0 r = 1 Lo w er b ounds obtained b y minimizing c T x ov er S r -1 -0.9420 Upp er b ounds obtained by minimizing c T x ov er I r -0.9105 -0.9420 T able 2: Our low er and upp er b ounds on the optimal v alue of the R-LD-LP in Example 2.12 . Figur e 8 b etter demonstr ates what the SDP is achieving at r = 1 . The set I 1 ( E 1 ) in Figur e 8a is the set of p oints in P that land in the set E 1 of Figur e 8b after one applic ation of G . Both E 1 and I 1 ( E 1 ) ar e by c onstruction invariant inner appr oximations to S . But as exp e cte d, E 1 ⊆ I 1 ( E 1 ) , which is why I 1 ( E 1 ) is the inner appr oximation of inter est at r = 1 . Note also that the el lipsoid E 1 that the SDP ﬁnds at r = 1 (Figur e 8b ) is very diﬀer ent fr om the el lipsoid E 0 than the SDP ﬁnds at r = 0 (Figur e 7a ). 3. Robust to uncertain and time-va rying linear dynamics linea r p rogramming. In the theory of robust con trol, there has lon g b een an in terest in analyzing the b ehavior of dynamical systems whose parameters are not exactly kno wn and can v ary in time. This is motiv ated by the fact that in man y practical applications, the physical or so cial dynamics of in terest are hard to mo del exactly and are sub ject to external disturbances that v ary with time. W e consider one of the most widely-studied linear mo dels that captures b oth parameter uncertain t y and dep endence on time (see, e.g., [ 28 ], [ 48 ] and references therein). In this mo del, one is given s real n × n matrices G 1 , . . . , G s and assumes that the true linear dynamics are giv en by a matrix in their con v ex hull conv { G 1 , . . . , G s } . (This is a p olytop e in the space of s × s matrices whose extreme p oin ts are giv en.) Moreov er, at each iteration, a diﬀerent matrix from this p olytop e can gov ern the dynamics. This leads to the following uncertain and 24 A. A. AHMADI AND O. G ¨ UNL ¨ UK time-v arying dynamical system (3.1) x k +1 ∈ conv { G 1 , . . . , G s } x k , where con v { G 1 , . . . , G s } x : = { P s i =1 λ i G i x | λ i ≥ 0 , P s i =1 λ i = 1 } . Note that this mo del en- compasses a nonlinear, time-inv ariant system x k +1 = g ( x k ) , where g : R n → R n satisﬁes (3.2) ∀ x ∈ Ω , g ( x ) ∈ conv { G 1 x, . . . , G s x } for the relev ant subset Ω of R n (see ( 1.1 )). Moreov er, this mo del captures uncertain t y in the nonlinear dynamics as well, as the dynamical system can b e gov erned b y any map g that satisﬁes ( 3.2 ). In this section, w e are interested in studying linear programs that must remain robust against suc h a dynamical sys tem. More precisely , a r obust to unc ertain and time-varying line ar dynamics line ar pr o gr am (R-UTVLD-LP) is an optimization problem of the form min x 0 ∈ R n  c T x 0 : x k ∈ P for k = 0 , 1 , 2 , . . . , u.t.d. x k +1 ∈ conv { G 1 , . . . , G s } x k  , (3.3) where P = { x ∈ R n | Ax ≤ b } is a giv en p olyhedron. The input to this problem is fully deﬁned b y A ∈ R m × n , b ∈ R m , c ∈ R n , G 1 , . . . , G s ∈ R n × n . It is not hard to see that an R-UTVLD-LP can b e equiv alently formulated as min x 0 ∈ R n  c T x 0 : x k ∈ P for k = 0 , 1 , 2 , . . . , u.t.d. x k +1 ∈ { G 1 x k , . . . , G s x k }  . Indeed, it is straightforw ard to chec k that for any in teger k ≥ 1, a p oin t x ∈ P leav es P by some pro duct of length k out of the matrices in conv { G 1 , . . . , G s } , if and only if it leav es P b y some pro duct of length k out of the matrices in { G 1 , . . . , G s } . Let G : = { G 1 , . . . , G s } and let G k denote the set of all s k matrix pro ducts of length k (with G 0 consisting only of the identit y matrix by con v en tion). Let G ∗ = ∪ ∞ k =0 G k b e the set of all ﬁnite pro ducts from G . An R-UTVLD-LP can then b e reformulated as the follo wing linear program with a countably inﬁnite n um b er of constraints: (3.4) min x ∈ R n { c T x : Gx ∈ P , ∀ G ∈ G ∗ } . Note that an R-LD-LP is a sp ecial case of an R-UTVLD-LP with s = 1. Throughout this section, we denote the feasible set of an R-UTVLD-LP b y (3.5) S := ∞ \ k =0 { x ∈ R n | AGx ≤ b, ∀ G ∈ G k } . ROBUST-TO-D YNAMICS OPTIMIZA TION 25 Clearly , the statemen t of Theorem 2.1 still applies to this set. Indeed, S is closed and con v ex as an inﬁnite in tersection of closed conv ex sets, and, b y deﬁnition, inv ariant under m ultiplication by G 1 , . . . , G s . Moreo ver, S is not alw ays p olyhedral even when s = 1, and testing mem b ership of a given p oin t to S is NP-hard already when s = 1. Our goal here will be to study tractable outer and inner approximations to S , and to extend some of the statemen ts we pro v ed for R-LD-LPs to this more intricate setting. 3.1. Outer appro ximations to S . Let (3.6) S r := r \ k =0 { x ∈ R n | AGx ≤ b, ∀ G ∈ G k } denote the set of p oin ts that remain in P under all matrix pro ducts of length up to r . It is clear that these sets provide p olyhedral outer approximations to S : S ⊆ . . . S r +1 ⊆ S r ⊆ . . . ⊆ S 2 ⊆ S 1 ⊆ S 0 = P . Hence, by solving LPs that minimize c T x o v er S r , w e obtain a nondecreasing and con v ergen t sequence of lo wer b ounds on the optimal v alue of an R-UTVLD-LP . W e lea v e it to the reader to chec k that the statement of Lemma 2.3 still holds with an almost iden tical pro of. This giv es us a wa y of c hec king ﬁnite termination of con v ergence of the sets S r to S . W e now need to generalize the notion of the sp ectral radius to several matrices. Deﬁnition 3.1 (Rota and Strang [ 45 ]). Given a set of n × n matric es G = { G 1 , . . . , G s } , their join t sp ectral radius (JSR) is deﬁne d as 12 ρ ( G ) := lim k →∞ max σ ∈{ 1 ,...,s } k || G σ 1 . . . G σ k || 1 /k . The JSR c haracterizes the maxim um growth rate that can b e obtained b y taking long pro ducts out of the matrices G 1 , . . . , G s in arbitrary order. Note that when s = 1, it coincides with the sp ectral radius. This can b e seen e.g., via the Gelfand’s formula for the sp ectral radius. W e observ e that the statemen ts of Prop ositions 2.4 , 2.5 , and 2.6 are still v alid (with ρ ( G ) replaced with ρ ( G )), as they even apply to the sp ecial case of an R-UTVLD-LP with s = 1. These prop ositions, together with the construction in the pro of of part (ii) of Theorem 2.1 , demonstrate that none of the three assumptions in the follo wing theorem can b e remov ed. Theo rem 3.2. L et G = { G 1 , . . . , G s } . If ρ ( G ) < 1 , P is b ounde d, and the origin is in the interior of P , then S = S r for some inte ger r ≥ 0 . Pr o of. Let ˆ ρ = ρ ( G )+1 2 < 1. It follows (see, e.g., [ 45 ], [ 11 , Lemma I I]) that there exists a norm f : R n → R suc h that for any α ≥ 0 and an y x ∈ R n , f ( x ) ≤ α ⇒ f ( G i x ) ≤ α · ˆ ρ , ∀ i ∈ { 1 , . . . , s } . As P contains the origin in its interior and is b ounded, there exists α 2 > α 1 > 0 suc h that { x ∈ R n | f ( x ) ≤ α 1 } ⊆ P ⊆ { x ∈ R n | f ( x ) ≤ α 2 } . Hence, any p oin t in P , once multiplied by an y matrix pro duct of length r = log( α 1 /α 2 ) log( ˆ ρ ) , lands in the set { x ∈ R n | f ( x ) ≤ α 1 } . As { x ∈ R n | f ( x ) ≤ α 1 } ⊆ S , the result follows. 12 The JSR is indep endent of the norm used in this deﬁnition. 26 A. A. AHMADI AND O. G ¨ UNL ¨ UK W e remark that our pro of ab o v e did not use the fact that P was a p olytope and w ould hold if P were instead any compact set. The reason this pro of was noticeably simpler than that of Theorem 2.7 is that w e did not analyze how large r can b e. W e did not do so b ecause of tw o reasons: (i) the sublevel sets of the norm f in the ab ov e pro of may not b e simple sets lik e ellipsoids that are amenable to algorithmic analysis, and (ii) ev en if r is small, the n umber of inequalities describing the set S r can b e as large as P r k =0 ms k , a quantit y which gro ws very quic kly when s ≥ 2. W e empirically observ e, how ev er, that the ﬁrst few lev els of this hierarch y often pro vide high-qualit y lo w er b ounds on the optimal v alue of an R-UTVLD-LP . W e can c hec k this b y computing upp er b ounds on the optimal v alue via a pro cedure that w e describ e in the next subsection. Theorem 3.2 as well as some of the theorems in the remainder of this section require the assumption that ρ ( G ) < 1. While algorithmic decidability of this condition is currently unkno wn [ 19 ], there is a large bo dy of literature on the computation of (arbitrarily tight) upp er b ounds on the JSR, whic h can b e utilized to verify this assumption; see e.g. [ 28 ], [ 42 ], [ 17 ], [ 3 ] and references therein. In fact, w e presen t a hierarch y of SDP-based suﬃcient conditions for c hec king this assumption in the next subsection (see Theorem 3.4 ), which happ ens to also b e useful for ﬁnding inner approximations to the feasible set of an R-UTVLD-LP . 3.2. Inner app ro ximations to S . In this subsection, w e generalize the results of Section 2.2 to the case of R-UTVLD-LPs. Recall our notation S from ( 3.5 ) for the feasible set of an R- UTVLD-LP , and let us k eep our notation P , G k , and S r from the previous subsection. Let E ⊆ P be an y con v ex set that con tains the origin in its interior and is in v ariant under m ultiplication by G 1 , . . . , G s . Since E is conv ex, it must also b e inv ariant under the dynamics in ( 3.1 ). Deﬁne S − 1 := R n and for an y integer r ≥ 0, let (3.7) I r ( E ) = S r − 1 ∩ { x ∈ R n | Gx ∈ E , ∀ G ∈ G r } . Note that I 0 ( E ) = E by deﬁnition. With this notation, the reader can verify that Lemma 2.8 , Lemma 2.9 , and Corollary 2.10 extend, with almost identical pro ofs, to the cas e where the single matrix G is replaced by the set of matrices G = { G 1 , . . . , G s } . W e summarize these results in the next lemma. Lemma 3.3. L et E ⊆ P b e c onvex 13 and invariant with r esp e ct to G = { G 1 , . . . , G s } . The sets I r ( E ) in ( 3.7 ) satisfy the fol lowing pr op erties: I r ( E ) ⊆ S , and I r ( E ) ⊆ I r +1 ( E ) for al l r ≥ 0 . Mor e over, if P is b ounde d, ρ ( G ) < 1 , and E c ontains the origin in its interior, then ther e exists a nonne gative inte ger t such that S = I r ( E ) for al l r ≥ t . In words, Lemma 3.3 states that the sets I r ( E ) pro vide an impro ving sequence of inner appro ximations to S and coincide with S in ﬁnite time. 13 W e ask that E b e conv ex, so its inv ariance with resp ect to G would imply its inv ariance with respect to the matrices in conv ( G ). It is easy to see that in general, if a set T is inv ariant under G , then conv ( T ) is inv ariant under conv ( G ). ROBUST-TO-D YNAMICS OPTIMIZA TION 27 3.2.1. Computation of I r ( E ) . The construction of the sets I r ( E ) requires access to a con v ex in v ariant set E ⊆ P . A non trivial c hallenge here is that unlike the case of a single matrix (Section 2.2.1 ), it is p ossible to hav e ρ ( G ) < 1 and y et not hav e an ellipsoid that is in v ariant under the action of the m atrices G 1 , . . . , G s . F or example, the matrices G 1 = γ  1 0 1 0  , G 2 = γ  0 1 0 − 1  ha v e JSR less than one for γ ∈ [0 , 1) , but only admit a common inv ariant ellipsoid for γ ∈ [0 , 1 √ 2 ] [ 6 ]. It turns out ho w ev er, that if the JSR is less than one, then there is al- w a ys an in v arian t set whic h is the in tersection of a ﬁnite num ber of ellipsoids. Moreo v er, these ellipsoids can b e found via semideﬁnite programming. Theo rem 3.4 (see Theorem 6.1 and Theorem 2.4 of [ 4 ]). L et G = { G 1 , . . . , G s } b e a set of n × n matric es. Then, for any inte ger l ≥ 1 , if ρ ( G ) ≤ 1 2 l √ n , ther e exist s l − 1 r e al symmetric matric es H π , wher e π ∈ { 1 , . . . , s } l − 1 is a multi-index, such that (3.8) H π ≻ 0 ∀ π ∈ { 1 , . . . , s } l − 1 , G T j H iσ G j ⪯ H σ j , ∀ σ ∈ { 1 , . . . , s } l − 2 , ∀ i, j ∈ { 1 , . . . , s } . Conversely, existenc e of a set of symmetric matric es H π that satisfy the semideﬁnite c on- str aints in ( 3.8 ) strictly 14 implies that ρ ( G ) < 1 . Mor e over, if ( 3.8 ) is satisﬁe d, then for any sc alar α ≥ 0 , the set (3.9) F α : =  x ∈ R n | x T H π x ≤ α, ∀ π ∈ { 1 , . . . , s } l − 1  is invariant under multiplic ation by G 1 , . . . , G s . R emark 3.5. By con v en tion, when l = 1, the decision v ariable in ( 3.8 ) is just a single matrix H and the constrain ts in ( 3.8 ) should read H ≻ 0 , G T j H G j ⪯ H , ∀ j ∈ { 1 , . . . , s } . In the case where l = 2, one should solve ( 3.8 ) with the conv ention that { 1 , . . . , s } 0 is the empt y set. This means that the decision v ariables are H 1 , . . . , H s and the constrain ts are H 1 ≻ 0 , . . . , H s ≻ 0 , G T j H i G j ⪯ H j , ∀ i, j ∈ { 1 , . . . , s } . Pr o of of The or em 3.4 . The pro of of this theorem app ears in [ 4 ], except for the part ab out in v ariance of the sets F α , which w e include here for completeness. W e need to show that the constrain ts in ( 3.8 ) imply x ∈ F α ⇒ G j x ∈ F α , ∀ j = 1 , . . . , s. 14 If ρ ( G ) < 1 2 l √ n , the constraints in ( 3.8 ) will indeed b e strictly feasible as one can apply the ﬁrst part of this theorem to β G : = { β G 1 , . . . , β G s } for β > 1 and small enough. (Note that the JSR is a con tinuous function of the entries of G [ 28 ] and satisﬁes the homogeneity relation ρ ( β G ) = β ρ ( G ).) 28 A. A. AHMADI AND O. G ¨ UNL ¨ UK Let ¯ x ∈ F α and deﬁne a function W : R n → R as W ( x ) := max π ∈{ 1 ,...,s } l − 1 { x T H π x } . By deﬁnition of F α , W ( ¯ x ) ≤ α . F urthermore, from the second set of inequalities in ( 3.8 ), it is easy to see that W ( G j x ) ≤ W ( x ) , ∀ j = 1 , . . . , s and x ∈ R n . Indeed, ( 3.8 ) implies that ∀ σ ∈ { 1 , . . . , s } l − 2 , ∀ i, j ∈ { 1 , . . . , s } and ∀ x ∈ R n , x T G T j H iσ G j x ≤ max ˆ σ ∈{ 1 ,...,s } l − 2 , ˆ j ∈{ 1 ,...,s } x T H ˆ σ ˆ j x = W ( x ) . W e hence deduce that W ( G j ¯ x ) ≤ W ( ¯ x ) ≤ α, for j = 1 , . . . , s, and so G j ¯ x ∈ F α for j = 1 , . . . , s. Going back to the computation of the conv ex in v ariant set E ⊆ P , which is needed for the construction of the inner appro ximations I r ( E ) in ( 3.7 ), w e ﬁrst ﬁnd the smallest in teger l ≥ 1 for which the SDP in ( 3.8 ) is feasible. (Note that we nev er need to compute the JSR.) Once this is done, for an y ﬁxed α ≥ 0, the set F α in ( 3.9 ) provides us with a conv ex and in v ariant set. W e no w need to ﬁnd a small enough ¯ α > 0 such that F ¯ α ⊆ P . A simple w a y of doing this is to require that one ellipsoid, sa y the ﬁrst, b e in the p olytop e. With this approach, ¯ α can b e computed b y follo wing the pro cedure describ ed in Step 2 of the pro of of Theorem 2.7 , whic h only requires matrix inv ersion. With ¯ α and the matrices { H π } ﬁxed, consider the following sequence of con v ex quadratic programs, minimize x ∈ R n c T x s.t. ( Gx ) T H π ( Gx ) ≤ ¯ α, ∀ G ∈ G r , ∀ π ∈ { 1 , . . . , s } l − 1 (3.10) AGx ≤ 1 , ∀ G ∈ r − 1 [ k =0 G k , indexed b y an integer r ≥ 0. The feasible sets of these optimization problems are the sets I r ( E ) as deﬁned in ( 3.7 ) with E = F ¯ α . As I r ( E ) ⊆ S for all r ≥ 0, the optimal v alues of these conv ex quadratic programs are upp er b ounds on the optimal v alue of the R-UTVLD- LP . Lemma 3.3 further implies that these upp er b ounds monotonically impro v e with r and reac h the optimal v alue of the R-UTVLD-LP in a ﬁnite num b er of steps. While this approach already achiev es ﬁnite conv ergence, there is muc h ro om for improv emen t as the in v ariant set E is ﬁxed throughout the iterations and is designed without taking in to consideration the ob jective function. 3.2.2. Computation of imp roved inner app ro ximations. Our goal no w is to ﬁnd in v arian t sets E r that result in the sets I r ( E r ) in ( 3.7 ) that b est approximate the feasible S of an R- UTVLD-LP in the direction of its ob jective function. T o do this, w e ﬁrst ﬁnd the smallest in teger l for which the SDP in ( 3.8 ) is feasible. W e ﬁx this num ber l throughout. Our sets E r , for r = 0 , 1 , . . . , will then b e given by E r =  z ∈ R n | z T H π ,r z ≤ 1 , ∀ π ∈ { 1 , . . . , s } l − 1  , (3.11) ROBUST-TO-D YNAMICS OPTIMIZA TION 29 where the symmetric matrices H π ,r are optimal solutions to the follo wing optimization prob- lem: minimize x ∈ R n ,H π ∈ S n × n c T x (3.12) s.t. H π ≻ 0 , ∀ π ∈ { 1 , . . . , s } l − 1 , G T j H iσ G j ⪯ H σ j , ∀ σ ∈ { 1 , . . . , s } l − 2 , ∀ i, j ∈ { 1 , . . . , s } , (3.13) ∀ z ∈ R n , z T H 1 ... 1 z ≤ 1 ⇒ Az ≤ 1 , (3.14) ( Gx ) T H π ( Gx ) ≤ 1 , ∀ π ∈ { 1 , . . . , s } l − 1 , ∀ G ∈ G r , (3.15) AGx ≤ 1 , ∀ G ∈ r − 1 [ k =0 G k . (3.16) Our Remark 3.5 regarding notation still applies here. Note that constraint ( 3.13 ) imp oses that the set E r in ( 3.11 ) b e inv ariant under the dynamics in ( 3.1 ). Constraint ( 3.14 ) forces one of the ellipsoids to b e within the p olytop e, which implies that the intersection E r of all ellipsoids will b e in the polytop e (this is obviously only a suﬃcient condition for E r ⊆ P ). W e remark here that c ho osing H 1 ... 1 to feature in this constraint is without loss of generality; as H 1 ... 1 is a v ariable of the problem, the optimization problem will naturally pick the “b est” ellipsoid to constrain to b e in the p olytop e. Constrain ts ( 3.15 ) and ( 3.16 ) force the p oint x to land in E r under all pro ducts of length r without lea ving P b efore time r . Once this optimization problem is solved to obtain an optimal solution x r and { H π ,r } , our inner approximation to S at step r will b e the set (3.17) I r ( E r ) = n z ∈ R n | ( Gz ) T H π ,r ( Gz ) ≤ 1 , ∀ π ∈ { 1 , . . . , s } l − 1 , ∀ G ∈ G r , AGz ≤ 1 , ∀ G ∈ r − 1 [ k =0 G k o , and x r will serve as our candidate sub optimal solution to R-UTVLD-LP . Just as we did in Section 2.2.2 , w e next sho w that b y a reparameterization, the ab o ve optimization problem can b e cast as an SDP . Theo rem 3.6. Supp ose ρ ( G ) < 1 and the set P = { x ∈ R n | Ax ≤ 1 } is b ounde d. L et a i denote the tr ansp ose of the i -th r ow of the matrix A ∈ R m × n and c onsider the fol lowing 30 A. A. AHMADI AND O. G ¨ UNL ¨ UK semideﬁnite pr o gr am: minimize x ∈ R n ,Q π ∈ S n × n c T x s.t. Q π ≻ 0 , ∀ π ∈ { 1 , . . . , s } l − 1 , (3.18) G j Q σ j G T j ⪯ Q iσ , ∀ σ ∈ { 1 , . . . , s } l − 2 , ∀ i, j ∈ { 1 , . . . , s } , a T j Q 1 ... 1 a j ≤ 1 , ∀ j ∈ { 1 , . . . , m }  Q π Gx ( Gx ) T 1  ⪰ 0 , ∀ G ∈ G r , ∀ π ∈ { 1 , . . . , s } l − 1 , AGx ≤ 1 , ∀ G ∈ r − 1 [ k =0 G k . Then, (i) the optimal values of pr oblems ( 3.12 ) and ( 3.18 ) ar e the same, the optimal ve ctors x r in the two pr oblems ar e the same, and the optimal matric es H π ,r and Q π ,r ar e r elate d via Q π ,r = H − 1 π ,r . (ii) the optimal values of the SDPs in ( 3.18 ) pr ovide upp er b ounds on the optimal value of the R-UTVLD-LP, ar e nonincr e asing with r , and r e ach the optimal value of the R-UTVLD-LP in a ﬁnite numb er of steps ¯ r . Mor e over, any optimal solution x r to the SDP in ( 3.18 ) with r ≥ ¯ r is an optimal solution to R-UTVLD-LP. Pr o of. The pro of of part (i) uses the same exact ideas as the pro of of part (i) of Theo- rem 2.11 (Sc h ur complemen ts and p olar duality of p olytop es and ellipsoids) and is left to the reader. In particular, this pro of w ould use Sc h ur complements to show that G T j H iσ G j ⪯ H σ j ⇐ ⇒ G j H − 1 σ j G T j ⪯ H − 1 iσ . W e now prov e part (ii). The statement that the optimal v alue of ( 3.18 ) is an upp er b ound on the optimal v alue of the R-LD-LP follo ws from the fact that in view of part (i), the last tw o sets of constrain ts of this SDP are constraining the optimal solution x r to b e in I r ( E r ), as deﬁned in ( 3.17 ). W e know that I r ( E r ) ⊆ S , ∀ r ≥ 0 as p oints in I r ( E ) land in the inv ariant set E r ⊆ P (cf. ( 3.11 )) in r steps without lea ving P b efore time r . T o see the claim ab out the monotonic improv ement of our upp er b ounds, observ e that if x, { Q π } are feasible to the SDP in ( 3.18 ) at level r , then they are also feasible to the SDP at lev el r + 1. This is b ecause w e hav e the inclusion E : =  z ∈ R n | z T Q − 1 π z ≤ 1 , ∀ π ∈ { 1 , . . . , s } l − 1  ⊆ P , b y the third set of constraints in ( 3.18 ), and the implication Gx ∈ E , ∀ G ∈ G r ⇒ Gx ∈ E , ∀ G ∈ G r +1 b y inv ariance of E under { G 1 , . . . , G s } as enforced by the second set of constraints in ( 3.18 ). W e now show that there exists an integer ¯ r ≥ 0 such that the optimal v alue c r to the SDP at level r is equal to the optimal v alue c ∗ of the R-UTVLD-LP for all r ≥ ¯ r . Let E = E 0 as deﬁned in ( 3.11 ). Observ e that the set E so deﬁned satisﬁes the assumptions of Lemma 3.3 ROBUST-TO-D YNAMICS OPTIMIZA TION 31 and hence there exists an integer ¯ r ≥ 0 suc h that S = I r ( E ) for all r ≥ ¯ r . Consequently , the optimal v alue of the conv ex quadratic program in ( 3.10 ) with H π = H π , 0 and ¯ α = 1 is equal to c ∗ for an y r ≥ ¯ r . As this optimal v alue is an upp erbound on the optimal v alue of the SDP in ( 3.18 ) for an y r ≥ 0 (indeed, H π , 0 is alwa ys feasible to ( 3.18 )), the claim follows. Finally , as an y optimal solution x r to the SDP at lev el r ≥ ¯ r satisﬁes c T x r = c ∗ and b elongs to S , it m ust b e an optimal solution to the R-UTVLD-LP as well. W e end with a n umerical example. Example 3.7. Consider an R-UTVLD-LP deﬁne d by the fol lowing data: A =       1 0 − 1 . 5 0 0 1 0 − 1 1 1       , b =       1 1 1 1 1       , c =  0 . 5 1  , G 1 = α  − 1 − 1 − 4 0  , and G 2 = α  3 3 − 2 1  , with α = 0 . 254 . F or this value of α (and in fact for any α ≥ 0 . 252 ), ther e is no el lipsoid that is invariant under multiplic ation by the p air G = { G 1 , G 2 } . This c an b e se en by observing that the SDP in ( 3.8 ) is infe asible when l = 1 . However, fe asibility of this SDP with l = 2 shows that ther e ar e two el lipsoids whose interse ction is invariant under the action of G 1 and G 2 . 15 In T able 3 , we give upp er and lower b ounds on the optimal value of this R-UTVLD-LP. T o obtain the lower b ounds, we minimize c T x over the sets S r in ( 3.6 ) for r = 0 , 1 , 2 . T o obtain the upp er b ounds, we solve the SDP in ( 3.18 ) for l = 2 and r = 0 , 1 , 2 . F or the c onvenienc e of the r e ader, we write out this SDP ( a T j her e denotes the j -th r ow of the matrix A ): minimize x ∈ R 2 ,Q 1 , 2 ∈ S 2 × 2 c T x (3.19) s.t. Q 1 ≻ 0 , Q 2 ≻ 0 , G 1 Q 1 G T 1 ⪯ Q 1 , G 2 Q 2 G T 2 ⪯ Q 1 , G 1 Q 1 G T 1 ⪯ Q 2 , G 2 Q 2 G T 2 ⪯ Q 2 ,  Q i Gx ( Gx ) T 1  ⪰ 0 , ∀ G ∈ G r , i = 1 , 2 , a T j Q 1 a j ≤ 1 , j = 1 , . . . , 5 , AGx ≤ 1 , ∀ G ∈ G k , k = 0 , . . . , r − 1 . F r om T able 3 , we note that as exp e cte d, our se quenc e of upp er b ounds (r esp. lower b ounds) ar e nonincr e asing (r esp. nonde cr e asing). Though we know that these b ounds must c onver ge to the optimal value of our R-UTVLD-LP in ﬁnite time, c onver genc e has not o c curr e d in this example in 3 iter ations. Inde e d, the gap b etwe en the upp er b ound and lower b ound for r = 2 is quite smal l but stil l nonzer o. 15 F or α ≥ 0 . 256, we hav e ρ ( G 1 , G 2 ) > 1 , and hence no compact full-dimensional set can b e inv arian t under the action of G 1 and G 2 . The fact that ρ ( G 1 , G 2 ) > 1 can b e seen b y observing that p ρ ( G 1 G 2 ) is a low er b ound on ρ ( G 1 , G 2 ) [ 28 ] and that p ρ ( G 1 G 2 ) = 1 . 0029 when α = 0 . 256. 32 A. A. AHMADI AND O. G ¨ UNL ¨ UK r = 0 r = 1 r=2 Lo w er b ounds obtained b y minimizing c T x ov er S r -1.3333 -0.9374 -0.8657 Upp er b ounds obtained by minimizing c T x ov er I r -0.7973 -0.8249 -0.8417 T able 3: Our low er and upp er b ounds on the optimal v alue of the R-UTVLD-LP in Exam- ple 3.7 . In Figur e 9 , we have plotte d the outer appr oximations S r to the set S in dark gr ay, and the inner appr oximations I r ( E r ) to the set S in light gr ay. T o b e mor e sp e ciﬁc, let Q 1 ,r , Q 2 ,r b e optimal matric es to the SDP in ( 3.19 ) at level r . The sets I r ( E r ) that ar e depicte d ar e deﬁne d as: I r ( E r ) = n z ∈ R n | ( Gz ) T Q − 1 1 ,r ( Gz ) ≤ 1 , ( Gz ) T Q − 1 2 ,r ( Gz ) ≤ 1 , ∀ G ∈ G r , AGz ≤ 1 , ∀ G ∈ r − 1 [ k =0 G k o . In e ach subﬁgur e, we have also plotte d the optimal solutions achieve d by minimizing c T x over the inner and outer appr oximations to S . Note that as r incr e ases, the set S gets sandwiche d b etwe en these two appr oximations mor e and mor e tightly. (a) r = 0 (b) r = 1 (c) r = 2 Figure 9: Three lev els of our inner and outer approximations to the feasible set of the R- UTVLD-LP in Example 3.7 . 3.3. A Sp ecial Case: P ermutation Matrices. In this section, w e consider a sp ecial case of R-UTVLD-LP when the input matrices gov erning the dynamical system in ( 3.1 ) are giv en b y n ! p erm utation matrices: P =  Π ∈ R n × n | Π is a p erm utation matrix  . Recall that a p ermutation matrix is a square binary matrix with each row and each column con taining exactly one nonzero en try . The joint sp ectral radius of the set of p ermutation matrices is exactly 1. ROBUST-TO-D YNAMICS OPTIMIZA TION 33 A useful property of P is that its elemen ts form a group under matrix m ultiplication. Consequen tly , we ha v e (3.20) Π 1 , Π 2 ∈ P = ⇒ Π 1 Π 2 ∈ P , and therefore all ﬁnite pro ducts of matrices from P belong to P . Moreo v er, the conv ex hull of P has a simple linear description [ 14 ]: (3.21) con v( P ) = n Π ∈ R n × n | X j ∈ N Π ij = 1 ∀ i ∈ N , X i ∈ N Π ij = 1 ∀ j ∈ N , Π ij ≥ 0 ∀ i, j ∈ N o where, N = { 1 , . . . , n } . When the linear dynamics of an R-UTVLD-LP is giv en by p erm uta- tion matrices, using the same notation as in ( 3.4 ) and the implication in ( 3.20 ), the feasible set of the R-UTVLD-LP can b e written as S : = ∞ \ k =0  x ∈ R n | A Π x ≤ b, ∀ Π ∈ P k  =  x ∈ R n | A Π x ≤ b, ∀ Π ∈ P  . W e next show that this feasible set has a p olynomial-size p olyhedral description in extended space (and hence R-UTVLD-LP can b e reform ulated as polynomial-size linear program in this case). Theo rem 3.8. Consider the R-UTVLD-LP in ( 3.4 ) with P = { x ∈ R n | Ax ≤ b } , wher e A ∈ R m × n , and G = P . Then, the fe asible set S of the R-UTVLD-LP c an b e written as S = n x ∈ R n | X i ∈ N u k i + X j ∈ N w k j ≤ b k ∀ k ∈ { 1 , . . . , m } , u k i + w k j ≥ A ki x j ∀ i, j ∈ N , k ∈ { 1 , . . . , m } o . Pr o of. First note that S = ∩ m k =1 S k , where S k =  x ∈ R n | a T k Π x ≤ b k , ∀ Π ∈ P  . Here, a T k denotes the k th ro w of A . F or a given p oint ¯ x ∈ R n , we ha v e ¯ x ∈ S k if and only if the optimal v alue z k to the linear program max  a T k Π ¯ x | Π ∈ P  = max  a T k Π ¯ x | Π ∈ con v( P )  do es not exceed b k . Using LP duality and ( 3.21 ), we hav e z k = max  X i,j ∈ N A ki ¯ x j Π ij | X j ∈ N Π ij = 1 ∀ i ∈ N , X i ∈ N Π ij = 1 ∀ j ∈ N , Π ij ≥ 0 ∀ i, j ∈ N  = min  X i ∈ N u k i + X j ∈ N w k j | u k i + w k j ≥ A ki ¯ x j ∀ i, j ∈ N  , where u k , w k ∈ R n are the dual v ariables asso ciated with the constraints in the maximization problem. As b oth problem are feasible, this implies that ¯ x ∈ S k if and only if the minimization problem ab ov e has a solution ¯ u k and ¯ w k with X i ∈ N ¯ u k i + X j ∈ N ¯ w k j ≤ b k . 34 A. A. AHMADI AND O. G ¨ UNL ¨ UK Consequen tly , we ha v e S k = n x ∈ R n | X i ∈ N u k i + X j ∈ N w k j ≤ b k , u k i + w k j ≥ A ki x j ∀ i, j ∈ N , u k , w k ∈ R n o and as S = ∩ m k =1 S k , the pro of is complete. W e remark that more generally , whenever one has a R-UTVLD-LP problem with dynamics matrices G in such a w a y that a p olynomial-size LP (or SDP) based description of con v( G ∗ ) is av ailable, then one can use the duality arguments ab o v e to reform ulate the problem as a p olynomial-size LP (or SDP). 4. F uture directions and a broader agenda: optimization with dynamical systems con- straints. In this pap er, we studied robust-to-dynamics optimization (RDO) problems where the optimization problem is an LP and the dynamics is gov erned either by a known, or an unkno wn and time-v arying linear system. Even in these t w o settings, a num ber of questions remain open. F or example, is the problem of testing membership of a point to the feasible set of an R-LD-LP decidable? W e hav e shown that this problem is NP-hard in general and p olynomial-time solv able when the sp ectral radius of the linear map is b ounded aw ay from one. It follows from [ 5 ] that the general problem is decidable when n ≤ 3. In higher dimen- sions, how ev er, decidability is unknown even in the sp ecial case where the input p olyhedron is a single halfspace. This sp ecial case is related to the so-called Skolem-Pisot problem; see, e.g., [ 30 , 29 ] and references therein. In the context of R-UTVLD-LPs, when ρ ( G ) < 1 and the p olytope P con taining the origin in its interior, can one analyze the num b er of steps needed for our inner and outer appro ximations to coincide as w e did for the case of R-LD-LPs? The algorithmic analysis of RDO problems in volving other t yp es of optimization problems and dynamical systems is also left for future research. W e note that in full generality , basic questions around RDO can b e undecidable. F or instance, it follows from [ 27 , Corollary 1] that when the dynamics is as in Section 3 with num b er of matrices s equal to 2 and the dimension n equal to 9, and when the input set Ω is the complemen t of a hyperplane, then testing mem b ership of a given p oin t to the feasible set S of the RDO is undecidable. Similarly , one can observe that the same question is undecidable when the dynamical systems is given by degree-2 p olynomial diﬀerential equations and the input set Ω is an op en halfspace; see [ 26 , Theorem 26], and also [ 41 ] for a similar result regarding dynamical systems arising from no- regret learning in games. Despite these negative results, we b eliev e that future research can iden tify additional sp ecial cases of RDO problems that admit tractable algorithms. More generally , w e b eliev e that optimization problems that incorp orate “dynamic al sys- tems (DS) c onstr aints” in addition to standard mathematical programming constrain ts can b e of interest to the optimization communit y at large. An optimization problem with DS constrain ts is a problem of the form: (4.1) minimize f ( x ) sub ject to x ∈ Ω ∩ Ω DS . Here, the set Ω is the feasible set of a standard mathematical program and is describ ed in an explicit functional form; i.e., Ω : = { x ∈ R n | h i ( x ) ≤ 0 , i = 1 , . . . , m } for some scalar v alued ROBUST-TO-D YNAMICS OPTIMIZA TION 35 functions h i . The constraint set Ω DS ho w ev er is deﬁned implicitly and alwa ys in relation to a dynamical system given in explicit form ˙ x = g ( x ) (in contin uous time) or x k +1 = g ( x k ) (in discrete time) . The set Ω DS corresp onds to p oints whose future tra jectory satisﬁes a pr esp e ciﬁe d desir e d pr op erty over time . The optimization problem ( 4.1 ) with diﬀerent DS constraints can lo ok lik e any of the following: Optimize f o ver the p oints in Ω whose future tra jectories under the dynamical system g • sta y in Ω for all time (invarianc e) , • asymptotically ﬂow to the origin (asymptotic stability) , • nev er intersect a given set Θ (c ol lision avoidanc e) , • reac h a giv en set Θ in ﬁnite time (r e achability) , etc. One can imagine these problems arising from mathematical mo dels in sev eral domains. F or example, in rob otics, the version with asymptotic stabilit y could corresp ond to establishing ho w high a humanoid rob ot can raise one of its legs b efore losing balance; see Figure 6 in [ 37 ]. The version with collision a v oidance could corresp ond to the problem of computing the highest initial sp eed at whic h an autonomous vehicle can p erform a turn without going oﬀ track. Figure 10 gives an example of a t wo-dimensional optimization problem with DS con- strain ts. Here, the ob jective function is f ( x ) = − x 2 , the set Ω (plotted in blue) is deﬁned b y a linear inequalit y and a con v ex quadratic inequality , and the dynamical system is a cubic diﬀeren tial equation, g ( x ) = ( − x 2 + 3 x 1 x 2 , x 1 − 1 2 x 2 1 x 2 ) T , whose resulting vector ﬁeld is plotted with little orange arrows. A DS constraint can b e an y of the four items listed ab ov e with Θ b eing the red triangle. Figure 10: An illustration of an optimization problem with dynamical systems (DS) con- strain ts. W e w ould cautiously argue that the optimization communit y is by and large not used to thinking ab out constrain ts that are deﬁned implicitly via a dynamical system. At the other end of the sp ectrum, while the study of inv ariant sets, regions of attraction, reac hable sets, 36 A. A. AHMADI AND O. G ¨ UNL ¨ UK etc. is very natural for control theorists, their fo cus often is not on settings where one needs to optimize o v er these sets sub ject to a melding with mathematical programming constraints. An interesting future researc h agenda would b e to initiate a systematic algorithmic study of optimization problems with DS constraints. By this, we mean a rigorous complexity study of problems of this type, via, e.g., either a polynomial-time algorithm, or an NP- hardness/undecidabilit y result, or an approximation algorithm. As can b e seen from the table b elo w, this pap er has only co v ered the tip of the iceb erg when it comes to such problems. Indeed, a class of optimization problems with DS constraints is deﬁned by taking one element of each of the three columns of this table. The starred entries corresp ond to cases to whic h this pap er has con tributed partial results. Opt. Problem “ f , Ω ” Type of Dynamical System “ g ” DS Constrain t “ Ω DS ” Linear program* Linear* In v ariance* Con v ex quadratic program Linear and uncertain/sto c hastic Inclusion in region of attraction Semideﬁnite program Linear, uncertain, and time-v arying* Collision a v oidance Robust linear program Nonlinear (e.g., p olynomial) Reac habilit y 0/1 integer program Nonlinear and time-v arying Orbital stability P olynomial program Discrete/con tinu ous/hybrid of b oth Sto chastic stabilit y . . . . . . . . . 5. Ackno wledgments. The authors are grateful to Georgina Hall for several insigh tful discussions (particularly around the con ten t of Prop osition 2.6 and Section 3.2.2 ), and to t w o anon ymous referees whose comments hav e grately improv ed our manuscript. REFERENCES [1] A. A. Ahmadi, A. Chaudhr y, V. Sindhw ani, and S. Tu , Safely le arning dynamic al systems . h ttps: //tin yurl.com/safe- learning- april- 2023 , 2023. [2] A. A. Ahmadi and O. G ¨ unl ¨ uk , R obust-to-dynamics line ar pr ogr amming , in Pro ceedings of the 54th IEEE Conference on Decision and Con trol, 2015, pp. 5915–5919. [3] A. A. Ahmadi and R. M. Jungers , L ower b ounds on c omplexity of Lyapunov functions for switche d line ar systems , Nonlinear Analysis: Hybrid Systems, 21 (2016), pp. 118–129. [4] A. A. Ahmadi, R. M. Jungers, P. A. P arrilo, and M. Roozbehani , Joint sp e ctr al r adius and p ath-c omplete gr aph Lyapunov functions , SIAM Journal on Control and Optimization, 52 (2014), pp. 687–717. [5] S. Almagor, J. Ouaknine, and J. Worrell , The semialgebr aic orbit pr oblem , vol. 126 of Leibniz In ternational Pro ceedings in Informatics, Schloss Dagstuhl, 2019, pp. 6:1–6:15. [6] T. Ando and M.-h. Shih , Simultane ous contr actibility , SIAM Journal on Matrix Analysis and Applica- tions, 19 (1998), pp. 487–498. [7] N. A thanasopoulos and R. M. Jungers , Computing the domain of attr action of switching systems subje ct to non-c onvex c onstr aints , in Pro ceedings of the 19th International Conference on Hybrid Systems: Computation and Control, A CM, 2016, pp. 41–50. [8] A. Bar vinok , A Course in Convexity , vol. 54, American Mathematical So c., 2002. [9] A. Ben-T al, L. El Ghaoui, and A. Nemiro vski , R obust Optimization , Princeton Universit y Press, 2009. [10] A. Ben-T al and A. Nemiro vski , R obust optimization–metho dology and applic ations , Mathematical Programming, 92 (2002), pp. 453–480. ROBUST-TO-D YNAMICS OPTIMIZA TION 37 [11] M. A. Berger and Y. W ang , Bounde d semigroups of matric es , Linear Algebra and its Applications, 166 (1992), pp. 21–27. [12] D. Ber tsimas, D. B. Bro wn, and C. Caramanis , The ory and applic ations of r obust optimization , SIAM Review, 53 (2011), pp. 464–501. [13] D. Ber tsimas, D. P achamano v a, and M. Sim , R obust line ar optimization under gener al norms , Op er- ations Research Letters, 32 (2004), pp. 510–516. [14] G. Birkhoff , T r es Observaciones Sobr e el Algebr a Line al , Univ. Nac. T ucuman, Ser. A, 5 (1946), pp. 147– 154. [15] G. Bitsoris , Positively invariant p olyhe dr al sets of discr ete-time line ar systems , International Journal of Con trol, 47 (1988), pp. 1713–1726. [16] F. Blanchini , Set invarianc e in c ontr ol , Automatica, 35 (1999), pp. 1747–1767. [17] V. D. Blondel and Y. Nestero v , Computational ly eﬃcient appr oximations of the joint sp e ctr al r adius , SIAM Journal on Matrix Analysis and Applications, 27 (2005), pp. 256–272. [18] V. D. Blondel and N. Por tier , The pr esenc e of a zer o in an inte ger line ar r ecurr ent se quenc e is NP-har d to de cide , Linear algebra and its Applications, 351 (2002), pp. 91–98. [19] V. D. Blondel and J. N. Tsitsiklis , A survey of c omputational c omplexity r esults in systems and c ontr ol , Automatica, 36 (2000), pp. 1249–1274. [20] S. Boyd and L. V andenberghe , Convex Optimization , Cambridge Universit y Press, 2004. [21] I. Cornfeld, S. Fomin, and Y. Sinai , Er go dic Theory , Grundlehren der mathematischen Wis- sensc haften, Springer-V erlag New Y ork, 1982. [22] M. Dahleh, M. A. Dahleh, and G. Verghese , L e ctur es on dynamic systems and c ontr ol , (2004). [23] H. Djelassi, A. Mitsos, and O. Stein , R e c ent advanc es in nonc onvex semi-inﬁnite pro gr amming: Applic ations and algorithms , EUR O Journal on Computational Optimization, 9 (2021), p. 100006. [24] C. E. Garcia, D. M. Prett, and M. Morari , Mo del pr e dictive contr ol: the ory and pr actic e—a survey , Automatica, 25 (1989), pp. 335–348. [25] E. G. Gilber t and K. T. T an , Line ar systems with state and c ontr ol c onstr aints: the the ory and applic ation of maximal output admissible sets , IEEE T ransactions on Automatic control, 36 (1991), pp. 1008–1020. [26] E. Hainr y , De cidability and unde cidability in dynamic al systems , Research rep ort: inria-00429965, 2009, pp. 1–27. [27] V. Hala v a and M. Hir vensalo , Impr ove d matrix p air unde cidability r esults , Acta Informatica, 44 (2007), pp. 191–205. [28] R. M. Jungers , The Joint Sp e ctr al Radius: The ory and Applic ations , vol. 385 of Lecture Notes in Con trol and Information Sciences, Springer, 2009. [29] T. Karimov, E. Kelmendi, J. Nieuwveld, J. Ouaknine, and J. W orrell , The p ower of p ositivity , 06 2023, pp. 1–11, https://doi.org/10.1109/LICS56636.2023.10175758 . [30] T. Karimov, E. Kelmendi, J. Ouaknine, and J. Worrell , What’s De cidable Ab out Discr ete Line ar Dynamic al Systems? , Springer Nature Switzerland, Cham, 2022, pp. 21–38. [31] E. C. Kerrigan and J. M. Maciejowski , Invariant sets for c onstr aine d nonline ar discr ete-time systems with applic ation to fe asibility in model pr e dictive c ontrol , in Decision and Control, 2000. Pro ceedings of the 39th IEEE Conference on, vol. 5, IEEE, 2000, pp. 4951–4956. [32] H. Khalil , Nonline ar Systems , Prentice Hall, 2002. Third edition. [33] M. Kord a, D. Henrion, and C. N. Jones , Convex c omputation of the maximum contr ol le d invariant set for p olynomial c ontr ol systems , SIAM Journal on Con trol and Optimization, 52 (2014), pp. 2944–2969. [34] M. K. Kozlo v, S. P. T arasov, and L. G. Khachiy an , The p olynomial solvability of convex quadr atic pr o gr amming , U.S.S.R. Comput. Math. Math. Phys., 20 (1980), pp. 223––228. [35] B. Lega t, P. T abuada, and R. M. Jungers , Computing c ontr ol led invariant sets for hybrid systems with applic ations to mo del-pr e dictive c ontr ol , Av ailable at arXiv:1802.04522, (2018). [36] M. L ´ opez and G. Still , Semi-inﬁnite pr o gr amming , Europ ean Journal of Operational Research, 180 (2007), pp. 491–518. [37] A. Majumdar, A. A. Ahmadi, and R. Tedrake , Contr ol and veriﬁc ation of high-dimensional systems with DSOS and SDSOS pr o gr amming , in 53rd IEEE Conference on Decision and Con trol, IEEE, 2014, pp. 394–401. [38] I. Mezi ´ c and S. Wiggins , A metho d for visualization of invariant sets of dynamic al systems b ase d on the 38 A. A. AHMADI AND O. G ¨ UNL ¨ UK er go dic p artition , Chaos: An Interdisciplinary Journal of Nonlinear Science, 9 (1999), pp. 213–218. [39] J. M. Mul vey, R. J. V anderbei, and S. A. Zenios , R obust optimization of lar ge-sc ale systems , Op er- ations research, 43 (1995), pp. 264–281. [40] I. Niven , Irr ational Numb ers , The Carus Mathematical Monographs, The Mathematical Asso ciation of America, 1956. [41] G. P. Andrade, R. Frongillo, and G. Piliouras , No-r egr et le arning in games is Turing c omplete , in Pro ceedings of the 24th A CM Conference on Economics and Computation, Association for Computing Mac hinery , 2023, p. 111. [42] P. A. P arrilo and A. Jadbabaie , Appr oximation of the joint sp e ctr al r adius using sum of squares , Linear Algebra and its Applications, 428 (2008), pp. 2385–2402. [43] R. Reemtsen and J. R ¨ uckmann , Semi-Inﬁnite Pr o gr amming , Springer-Science, New Y ork, 1998. [44] P. Rost alski and B. Sturmfels , Dualities in c onvex algebr aic ge ometry , Av ailable at (2010). [45] G. C. Rot a and W. G. Strang , A note on the joint sp e ctr al r adius , Indag. Math., 22 (1960), pp. 379– 381. [46] J. O. R oyset , Set-c onver genc e and its applic ation: A tutorial , Set-V alued and V ariational Analysis, 28 (2020), pp. 707–732. [47] A. Schrijver , The ory of Line ar and Inte ger Pr ogr amming , John Wiley & sons, 1998. [48] R. Shor ten, F. Wir th, O. Mason, K. Wulff, and C. King , Stability criteria for switche d and hybrid systems , SIAM review, 49 (2007), pp. 545–592. [49] M. Sipser , Intr oduction to the The ory of Computation , v ol. 2, Thomson Course T ec hnology , 2006.

Robust-to-Dynamics Optimization

Original Paper

Comments & Academic Discussion

Leave a Comment

Original Paper

Related Papers

Comments & Academic Discussion

Leave a Comment