Model Predictive Tracking Control for Invariant Systems on Matrix Lie Groups via Stable Embedding into Euclidean Spaces

1 Model Predicti v e T racking Control for In v ariant Systems on Matrix Lie Groups via Stable Embedding into Euclidean Spaces Dong Eui Chang, Karmvir Singh Phogat, and Jongeun Choi, Senior Member , I EEE Abstract —Fo r controller design for systems on manifolds embedded in Euclidean space, it is con venient to utilize a theory that requires a si ngle global coordinate system on the ambient Euclidean space rather than multiple local charts on the manif old or coordinate-fr ee tools from differ ential geometry . In this article, we apply such a theory to design model predictiv e tracking controllers f or systems whose dynamics evo lve on manifolds and illustrate its efficacy with th e fully actuated rigid body attitude control system. Index T erms —Model predictiv e control, Matrix Lie gr oups, T racking control, Attit u de c ontrol. I . I N T R O D U C T I O N M ODEL predictive control (MPC), which requires solv- ing a constrained finite time-horiz o n optimal co ntrol problem , has been initially utilized mostly in slow proc e ss industries [1]. In contrast to con ventional co ntrol sche mes, MPC is prevalent in safety critical systems due to its ability to hand le state and con trol constraints fo r large-scale sys- tems [2], [3]. Due to the rise in computatio n al power and sophisticated algorithm s, several successful MPC imp lemen- tations have recently been rep o rted in various applicatio ns with fast dyna m ics, inclu d ing autonom o us vehicles [4], [5] and power electronics [6]. O bvio usly , MPC design ing strate- gies for co ntinuou s-time systems req uire lin earized discrete- time sy stems, accou n ting for the system dynam ics. Beca u se linearization an d discretization o f th e system dynamics ar e relativ ely daun tin g tasks o n m anifolds as compared to Eu- clidean spa ces, designin g MPC on ma n ifolds is a non trivial matter . First, a manifold cannot be entirely covered by o ne local coordin ate chart unless it is d iffeomorphic to a Euclidean space. As a resu lt, one needs to carry out coordinate changes when the system trajector y tr av erses throu gh multiple char ts. Second, linearizatio n and discretization of system d ynamics are both local app roxima tio ns, so th ese pr ocedur es re q uire use of local char ts as well. In gen eral, coordin ate chan g e is an expensiv e o peration in terms of computatio n tim e, and it m ay Dong Eui Chang an d Karmvir Singh Phoga t are with th e School of Electric al Engineeri ng, KAIST , Da ejeon, 34141, Korea . dechang@kaist. ac.kr, karmvir.p@gmail.com Jongeun Choi is with the School of Mechani cal Engineeri ng, Y onsei Uni versi ty , Seoul 03722, Kore a. jo ngeunchoi@yonse i.ac.kr This re search h as been in part supported by KAIST under grant s G04170001 a nd N11180231 and by the Mid-caree r Research Program through the National Researc h Foundation of Korea funded by the Ministry of Science and ICT under grant NRF-2018R1A2B600 8063. (Corresponding Author: K.S. Phogat, Co-correspon ding Author: J. Choi). (T o be) Published in IEEE Tra nsactions on Automatic Control. doi: 10.1109/T AC .2019.2946231 introdu c e fairly large discontin u ities to th e dyn amics du e to switchings of local cost functions for M PC that are defined chartwise on the manifold . As m any m echanical sy stem s such as aerial vehicles and rob otic systems ev olve on manif olds, it is e ssen tial to have a theor y that does not require switching of charts or u nconventional tools from geo metric contr o l theo ry . In this article, we propose a mo del pr edictive con troller design for systems o n m a n ifolds by employin g a stable em- bedding technique , which is summarized as follows: Given a control system on a manifold M , fir st emb e d the manifold into the Euclidean space R n and extend the system dy namics stably to the ambient Euclidean space, i.e., the system dyn amics are extended on R n in such a way tha t the manifo ld M becom es an in v ariant attracto r o f the mod ified dyn amics defined on R n . Since the system dynamics are now de fin ed in R n , we can ca r ry out linearizatio n in one single glo bal Eu clidean coordin ate set for R n and then discretization of th e linear ize d system in R n . Th e stable embedd in g techniqu e in creases th e dimension o f the system, and therefor e MPC for the ex- tended dynamics may be computationally mor e expensive than the dynamics in minimal co ordinates; howe ver , it simplifies linearization and discretization of th e system dynam ics to a large extent. This approach was su ccessfully applied for linear stabilizing /tracking contro ller design [7] an d structur e- preserving nu merical integration [8]. Recent attemp ts on MPC o n manifo lds may be f o und in [9], [10], which req u ire implicit representation of th e system dynamics or explicit constraints in the op timization to preserve the manifold structur e of SO(3) . In additio n , these schemes require switching of charts as the local contr o l law , which is neede d for stability , is defined in ch arts. In co ntrast to these conventional schemes, we take the aforemention ed stable embedd in g approach to design MPC for systems d efined on manifold s. In ou r study , we co nsider a class o f systems defined by fully actu ated left inv ariant vector field s on matr ix L ie group s and stably em bed the system dyn amics in Euclidean spaces. Subsequ ently , to track a refe r ence trajectory , time- varying tr acking er r or dyn amics are d efined in the amb ient Euclidean spa c e; th ose are linearized along the reference and simp lified u sing sym metry inv ariance, b o th in on e sing le global co ordina te system on the a m bient Eu clidean spa c e, and the lo cal stabilizability of the origin a l no nlinear trackin g error dynamics to zero is then readily established. For ap plying MPC, the error dyna mics are line a rized and discretized , and the stabilizability of th e discrete- time linear error d ynamics is also proven, for which a f undam ental sufficient con di- 2 tion is der iv ed in an inequality form that inv olves the two parameters: the d iscretization time step and the transversal stability p arameter that is intr o duced in the proc ess of stable embedd in g. Later, an M PC law is designed for the discrete- time linear error dynamics and is ap p lied to th e original nonlinear system. It is worth mentionin g at this junctu re that, to the b est of the auth ors’ knowledge, the issue of establishin g exponential stability of the time-varying system dynamics under the synthesized MPC control law rem ains o pen. Some results on the stability of the samp led data sy stems may be found in [11], [1 2]. T o demonstrate the efficac y of the propo sed MPC tech nique, we desig n a tr a cking controller for spacecraft attitude contro l dyn amics and c onduc t n umerical studies for various real-time scen arios, suc h as referen ce track- ing un der tig ht con trol con straints and noisy m e asurements, to illustrate the validity of the designed co ntroller . Nu merical simulations show that the MPC tracking contr o ller de signed using the discrete- tim e linear system in the Euclid ean space is robust to un modeled disturbance s and pr ovides good trackin g perfor mance when applied to the actu al non linear system. The stru cture of this article is as follows: W e establish stabilizability of the tracking erro r dyn amics fo r left inv ariant systems on matrix Lie gro ups and discuss the design proced ure for the MPC contro ller in Section II . Section III is dev oted to design of an MPC tracking con tr oller fo r a r igid body attitude control system. Numerical studies of the d esigned track ing controller f or attitud e dyn amics are in Section IV, fo llowed by o u r con clusions in Section V . I I . M P C O N M A T R I X L I E G RO U P S Let G be a m atrix L ie gro up of d imension m with I a s the group identity and g be the Lie alg ebra o f the Lie group G . Suppose that a controlled system dyn amics on th e matrix L ie group is given by a left in variant vector field: Σ : ( ˙ g = g ξ , ˙ ξ = f ( ξ , u ) , (1) where g ∈ G, ξ ∈ g and u ∈ R m . I t is assum ed that ∂ f ∂ u ( ξ , u ) is an in vertible lin ear map from R m to g for each ( ξ , u ) ∈ g × R m . Su ppose that the matrix L ie g roup G is embedded into a E uclidean space R n × n . The vector space R n × n is split into two orthog o nal subspaces g and g ⊥ such that R n × n = g ⊕ g ⊥ , where g ⊥ is the orth ogona l co mplement of g in R n × n under the E uclidean metric defined by h A, B i = trace( A T B ) fo r all A, B ∈ R n × n . In the subsequ ent discussion , we ref er to g as the parallel dir ection an d g ⊥ as the tran sversal direction, an d we defin e orth ogon a l projectio n maps fr om the amb ient space R n × n to g an d g ⊥ as R n × n ∋ v 7→ v k ∈ g , R n × n ∋ v 7→ v ⊥ ∈ g ⊥ . A d etailed discussion on left inv ariant systems and the differ- ential g eometric tools e m ployed in this article may be fou n d in [13]–[ 15]. Let us tur n to stably embed the system dyn a m ics (1) in to R n × n considerin g the following assumptio n : Assumption 1 . There exis ts a C 2 function R n × n ∋ x 7→ V ( x ) ≥ 0 ∈ R with the following pr operties: (A-i) V − 1 (0) = G . (A-ii) V ( xg ) = V ( x ) for all x ∈ R n × n , g ∈ G. (A-iii) ∇ 2 V ( I ) is po sitive definite in the transver sal direction, i.e., ∇ 2 V ( I ) · ( y , y ) > 0 for all y ∈ g ⊥ \{ 0 } . Since V attains its minimum value 0 at each poin t in G , ∇ V vanishes on G . Hen ce, th e system d y namics (1) can be extended to th e ambien t Euclid e an space R n × n × g as ˜ Σ : ( ˙ x = xξ − α ∇ V ( x ) , ˙ ξ = f ( ξ , u ) , (2) where α > 0 , x ∈ R n × n , ξ ∈ g and u ∈ R m . Let R ∋ t 7→ ( g 0 ( t ) , ξ 0 ( t )) ∈ G × g (3) be a re f erence state trajecto ry and R ∋ t 7→ u 0 ( t ) ∈ R m be the correspo nding contr ol trajec to ry of the system dynamics (1). Assumption 2. There exist consta nts β g max ≥ β g min > 0 such that the refer en ce trajectory R ∋ t 7→ g 0 ( t ) ∈ G satisfies β g min I  g 0 ( t )  g 0 ( t )  ⊤  β g max I for a ll t. The track ing erro r trajectory , d efined b y R ∋ t 7→ ( E ( t ) , Ξ( t )) : =  x ( t ) g − 1 0 ( t ) − I , ξ ( t ) − ξ 0 ( t )  ∈ R n × n × g such that ( E ( t ) , Ξ( t )) = 0 fo r all t ensur es, that the sy stem dynamics (2) is tracking the referen ce trajectory , i.e., x ( t ) = g 0 ( t ) , ξ ( t ) = ξ 0 ( t ) for all t . Theref ore, the refer ence tracking problem is translated to stab ilization of the error dynam ics to zero. The error dyn amics fo r a given ref erence trajecto ry is giv en as δ Σ : ( ˙ E =  g 0 + E g 0  Ξ g − 1 0 − α ∇ V ( g 0 + E g 0 ) g − 1 0 , ˙ Ξ = f (Ξ + ξ 0 , u ) − f ( ξ 0 , u 0 ) . (4) T o design a linear MPC c o ntroller, let us linearize the erro r dynamics (4) aroun d zero . The linearized erro r dyna m ics around zero is given by δ Σ ℓ : ( ˙ E = g 0 Ξ g − 1 0 − α  ∇ 2 V ( g 0 ) · ( E g 0 )  g − 1 0 , ˙ Ξ = ∂ f ∂ ξ ( ξ 0 , u 0 )Ξ + ∂ f ∂ u ( ξ 0 , u 0 ) δ u, (5) where δ u : = u − u 0 . Befor e simplifyin g the error dynam ic s ˙ E in (5) and splitting it further along the parallel and the transversal direction to gain mo re geometr ic in sight, let u s discuss som e key prop erties associated with the function V : Lemma 1. Under Assump tion 1, the fo llowing h o ld: (a) ∇ 2 V ( I ) · v k = 0 for all v k ∈ g . (b) ∇ 2 V ( I ) · v ⊥ ∈ g ⊥ for a ll v ⊥ ∈ g ⊥ . (c) ∇ 2 V ( g ) · ( v g ) =  ∇ 2 V ( I ) · v  g − 1  ⊤ for all v ∈ R n × n and g ∈ G . 3 Pr oof. 1) Note that ∇ V ( g ) = 0 for all g ∈ G. Ther efore, ∇ 2 V ( I ) · v k = d ds     s =0 ∇ V (exp ( sv k )) = 0 . 2) For any v ⊥ ∈ g ⊥ , ∇ 2 V ( I ) · v ⊥ ∈ g ⊥ if and only if  v k , ∇ 2 V ( I ) · v ⊥  = 0 fo r all v k ∈ g . Theref o re, using the fact that ∇ 2 V ( I ) is symmetric and then ap plying Lemma 1(a) leads to the fo llowing: D v k , ∇ 2 V ( I ) · v ⊥ E = D ∇ 2 V ( I ) · v k , v ⊥ E = 0 for all v k ∈ g . 3) For ar b itrary vectors w , v ∈ R n × n ,  w, ∇ 2 V ( g ) · ( v g )  = d dt d ds     t = s =0 V ( g + tv g + sw ) = d dt d ds     t = s =0 V ( I + tv + sw g − 1 ) =  wg − 1 , ∇ 2 V ( I ) · v  =  w,  ∇ 2 V ( I ) · v  ( g − 1 ) ⊤  . Therefo re, we con clude that ∇ 2 V ( g ) · ( v g ) =  ∇ 2 V ( I ) · v  ( g − 1 ) ⊤ . Employing th e prop erties discussed in Lemm a 1, it is straightfor ward to show that the linearized er ror dyna mics in (5) is transfo rmed to the f ollowing: ˙ E ⊥ = − α   ∇ 2 V ( I ) · E ⊥  g 0 g ⊤ 0  − 1  ⊥ , (6a) ˙ E k = g 0 Ξ g − 1 0 − α   ∇ 2 V ( I ) · E ⊥  g 0 g ⊤ 0  − 1  k , (6b) ˙ Ξ = ∂ f ∂ ξ ( ξ 0 , u 0 )Ξ + ∂ f ∂ u ( ξ 0 , u 0 ) δ u, (6c) where the linear er ror E in (5) has simplified a nd sp lit into the tr ansversal direction error R ∋ t 7→ E ⊥ ( t ) ∈ g ⊥ and the parallel d ir ection erro r R ∋ t 7→ E k ( t ) ∈ g . Theorem 1 . Und e r Assump tions 1 an d 2, the system dynamics (6a) is exponentially stable to zer o. Pr oof. Let us d efine a cand id ate L yapunov fun ction g ⊥ ∋ η 7→ V ( η ) : = 1 2  ∇ 2 V ( I ) · η , η  ∈ R . Then, using th e p roperties of V discussed in Lemma 1, we obtain: d V dt  E ⊥  =  ∂ V ∂ η ( E ⊥ ) , ˙ E ⊥  = − α  ∇ 2 V ( I ) · E ⊥ ,   ∇ 2 V ( I ) · E ⊥  g 0 g ⊤ 0  − 1  ⊥  = − α D ∇ 2 V ( I ) · E ⊥ ,  ∇ 2 V ( I ) · E ⊥  g 0 g ⊤ 0  − 1 E ≤ − α β g max  ∇ 2 V ( I ) · E ⊥ , ∇ 2 V ( I ) · E ⊥  ≤ − 2 αλ min β g max V  E ⊥  , where λ min is the sm allest eigenv alue o f the op e rator ∇ 2 V ( I ) restricted to g ⊥ . T h erefor e, V  E ⊥ ( t )  ≤ exp  − 2 αλ min β g max t  V  E ⊥ (0)  , and b y the definition of V , we know that λ min k E ⊥ k 2 ≤ 2 V ( E ⊥ ) ≤ λ max k E ⊥ k 2 , where λ max is the largest eig en value of the oper ator ∇ 2 V ( I ) restricted to g ⊥ . So , k E ⊥ ( t ) k ≤ r λ max λ min exp  − αλ max β g max t  k E ⊥ (0) k . This proves the assertion. Remark 1. The pr oof of Theor em 1 uses the bo unded ness of g 0 ( t )  g 0 ( t )  ⊤ fr om above, see Assump tion 2. The lin e a rized er ror dynamics ( 6) can be exponentially stabilized to zero by feedback. Consequently , the or ig inal non- linear tra c k ing e r ror dy n amics (4) is exponentially stabilizable. Theorem 2. Suppose Assumptions 1 and 2 h o ld. Then, for any two matrices K p , K d ∈ R n × n such that the ma trix  0 I K p K d  (7) is Hurwitz, the PD-like co ntr oller δ u =  ∂ f ∂ u ( ξ 0 , u 0 )  − 1 n [Ξ , ξ 0 ] + Y − ∂ f ∂ ξ ( ξ 0 , u 0 )Ξ o (8) with Y : = g − 1 0  − ˙ W + K p E k + K d ( g 0 Ξ g − 1 0 + W )  g 0 (9) and W : = − α   ∇ 2 V ( I ) · E ⊥  g 0 g ⊤ 0  − 1  k , stabilizes the con tr olled dyn a mics ( 6) exponen tia lly to zer o, wher e ˙ W in (9) can be e xpr essed in terms of state v ariables using (6a) a n d ˙ g 0 = g 0 ξ 0 . F u rthermor e, the contr oller u = u 0 + δ u exponentially stab ilizes (4) to zer o. Pr oof. Since the exp onential stability of the tr a nsversal dy - namics (6a) has been proved in Theorem 1 , the expon ential stability of the su b systems (6 b) and (6c) with the contro l law (8) re m ains to be p roved. Apply ing the contr oller (8) to subsystems (6 c) transf o rms the subsystems (6b) and (6c) to the following system of differential e q uations ˙ E k ˙ ˜ e ! =  0 I K p K d   E k ˜ e  , (10) where ˜ e : = g 0 Ξ g − 1 0 − α   ∇ 2 V ( I ) · E ⊥  g 0 g ⊤ 0  − 1  k . (11) Therefo re, the linear system ( 10) is expon entially stable to zero if the matrix (7) is Hurwitz stable. Since E ⊥ and ˜ e are exponentially stable, it follows fr o m ( 1 1), Theorem 1 and Assumption 2 that Ξ is a lso expon e ntially stable. By the L ya- punov linearization metho d, th is contro l law also exponen tially stabilizes (4) to zero. T his p roves the assertion. 4 T o apply th e MPC, let us discr e tize the linearized error dynamics (6) using Euler’ s method as E ⊥ k +1 = E ⊥ k − hα   ∇ 2 V ( I ) · E ⊥ k  g 0 ,k g ⊤ 0 ,k  − 1  ⊥ , (12a) E k k +1 = E k k + hg 0 ,k Ξ k g − 1 0 ,k − αh   ∇ 2 V ( I ) · E ⊥ k  g 0 ,k g ⊤ 0 ,k  − 1  k , (12b) Ξ k +1 = Ξ k + h ∂ f ∂ ξ ( ξ 0 ,k , u 0 ,k )Ξ k + h ∂ f ∂ u ( ξ 0 ,k , u 0 ,k ) δ u k , (12c) where h is a d iscretization step, and for a fun c tion R ∋ t 7→ Γ( t ) ∈ R n × n , Γ k : = Γ( k h ) . The stability of the MPC for d iscr e te-time systems re quires stabilizability of th ese discrete-time systems [16]. Theref ore, it is crucia l to p rove stabilizability of the discrete-time dynam ic s (12). First, we prove the exponentially stability o f the subsystem (12a) for an ap propr iate ch oice of th e discr e tization step h an d the parameter α . T hen, we pr oceed to the gener al c a se and establish stabilizab ility of the discrete-time dynamics ( 12). In a d dition, T heorem 3 establishes a re la tio n b etween the stabilizing param eter α and the d iscretization step leng th h that is cruc ial in implemen tation o f MPC. Theorem 3. Supp ose that Assumptions 1 and 2 hold . Then, the system dyn a mics ( 12a) is exponentially stable to zer o if the following h olds: 0 < αh < 2 λ min β 2 g min λ 2 max β g max , wher e λ min and λ max ar e the min imum and th e maximum eigen- values of the operator ∇ 2 V ( I ) restri cted to g ⊥ , respectively . Pr oof. Let us d efine a cand id ate L yapunov fun ction g ⊥ ∋ η 7→ V d ( η ) : =  ∇ 2 V ( I ) · η , η  ∈ R . Then, using the prop erties of V fro m Lemm a 1 g iv es V d  E ⊥ k +1  =  ∇ 2 V ( I ) · E ⊥ k +1 , E ⊥ k +1  =  ∇ 2 V ( I ) · E ⊥ k , E ⊥ k  − 2 αh D ∇ 2 V ( I ) · E ⊥ k ,  ∇ 2 V ( I ) · E ⊥ k  g 0 ,k g ⊤ 0 ,k  − 1 E + α 2 h 2     ∇ 2 V ( I ) · E ⊥ k  g 0 ,k g ⊤ 0 ,k  − 1    2 ∇ 2 V ( I ) ≤  1 − 2 αh λ min β g max + α 2 h 2 λ 2 max β 2 g min  V d  E ⊥ k  ≤   αh λ max β g min − 1  2 + 2 αh  λ max β g min − λ min β g max   V d  E ⊥ k  . Therefo re, the system dynam ics (1 2a) is expone ntially stable if 0 ≤  αh λ max β g min − 1  2 + 2 αh  λ max β g min − λ min β g max  < 1 which lead s to the following condition : 0 < αh < 2 λ min β 2 g min λ 2 max β g max . This p roves the assertion. Theorem 4. Supp ose that Assumptions 1 and 2 hold . Then, for a ny two matrices K p , K d ∈ R n × n such that the ma trix  I hI K p K d  (13) is Schur stable, the contr oller δ u k = 1 h  ∂ f ∂ u ( ξ 0 ,k , u 0 ,k )  − 1 n Y k −  I + h ∂ f ∂ ξ ( ξ 0 ,k , u 0 ,k )  Ξ k o (14) wher e Y k : = g − 1 0 ,k +1  K p E k k + K d  g 0 ,k Ξ k g − 1 0 ,k + W k  − W k +1  g 0 ,k +1 and W k : = − α   ∇ 2 V ( I ) · E ⊥ k  g 0 ,k g ⊤ 0 ,k  − 1  k , stabilizes the contr olled dyn amics (1 2) exponentially to zer o. Pr oof. Since the exp onential stability of the tr a nsversal dy - namics (12a) has b een proved in Th e orem 3, the exponen tial stability of the subsystems (12b) and ( 12c) with the control law (14) re m ains to b e proved. Applying the controller (1 4) to subsystems (12c) tran sforms the subsystems (1 2b) an d (12c) to the following system of d ifference equ ations  E k k +1 ˜ e k +1  =  I hI K p K d   E k k ˜ e k  (15) where ˜ e k : = g 0 ,k Ξ k g − 1 0 ,k − α   ∇ 2 V ( I ) · E ⊥ k  g 0 ,k g ⊤ 0 ,k  − 1  k . (16 ) Therefo re, the linear system ( 15) is expon entially stable to zero if the matrix (13) is Schur stable. Since E k k and ˜ e k are exponentially stable, it follows fr o m ( 1 6), Theorem 3 and Assumption 2 that Ξ is also expone ntially stable. This proves the assertion. Equipp e d with Theorem 4, we are in a position to design an MPC con trol law for the system dyn amics (12). MPC computes a static state feedba ck control law at each time instant by solvin g a co n strained finite horizo n discrete-time optimal co ntrol pro blem. A typical optimal co ntrol pro blem for a ho rizon N is to minimize a p e rforma n ce objective J ( E 0: N , Ξ 0: N , δ u 0: N − 1 ) : = N − 1 X k =0 ( k E k k 2 Q E + k Ξ k k 2 Q Ξ ) + N − 1 X k =0 k δ u k k 2 Q δu + k E N k 2 Q f E + k Ξ N k 2 Q f Ξ , (17) where Q E , Q Ξ , Q δu , Q f E , Q f Ξ ∈ R n × n are po siti ve semidefi- nite m atrices, while satisfying th e system dynamics (12) and the state and control constraints E k ∈ X E k for all k = 0 , . . . , N , Ξ k ∈ X Ξ k for all k = 0 , . . . , N , δ u k ∈ U k for all k = 0 , . . . , N − 1 , (18) where X E k , X Ξ k are the admissible state sets and U k is an admissible action set at eac h time instant k . 5 Concisely , an op tim al con trol δ u j | j at th e time in stant j for a fixed given state  E j | j , Ξ j | j  is ob tained by solving the following c onstrained d iscrete-time optima l contro l p roblem: minimize { δu j + i | j } N − 1 i =0 J  E j : j + N | j , Ξ j : j + N | j , δ u j : j + N | j  subject to      dynamics (12) f or k = j | j, . . . , j + N − 1 | j, constraints ( 18) for k = j | j, . . . , j + N − 1 | j,  E j | j , Ξ j | j  is fixed . (19) Then, the con trol law u ( t ) = u 0 ( t ) + δ u j | j is app lied to the sy stem (1) fo r t ∈ [ j h, ( j + 1) h [ , where j = 0 , 1 , 2 , · · · . Remark 2 . Note that the system d y namics (1 2) is expo- nentially sta b ilizable. Ther efor e, we design an exponentially stabilizing MPC law for the dyn amics (1 2) in Euclidean spaces that in turn stabilizes Euler’ s appr oximation of the err or dyn amics (4) exponen tially to zer o [17]. However , to the be st of the a uthors’ knowledge, the issue of establishing exponential stability of the time-varying sampled data system (4) under the synthesized MPC contr ol law r emains open. Some r esults on the stability o f the samp le d data systems ma y be fou n d in [11], [12]. I I I . A N I L L U S T R A T I V E E X A M P L E : T H E R I G I D B O DY C O N T RO L S Y S T E M Let us consider an example o f r igid body attitude d ynamics to d iscuss the the o ry d ev eloped in Section II. The rigid body attitude co ntrol system is defined b y ˙ R = R ˆ Ω , (20a) ˙ Ω = I − 1 ( I Ω × Ω) + I − 1 u, (20b) where R ∈ SO(3) (th e set of 3 × 3 real or thogon al matrices with determ inant 1) is a ro tation matrix that determ ines the attitude of the rig id b ody , Ω ∈ R 3 defines the an gular veloc ity of the r igid body ; u ∈ R 3 is the co ntrol torque; I is the momen t of inertial matrix of the rigid b o dy; an d the hat map ∧ maps R 3 vectors to 3 × 3 real ske w symmetr ic ma tr ices such that ˆ x y = x × y for all x, y ∈ R 3 with × as th e vector prod uct on R 3 . Note that the man ifold SO (3) × R 3 ⊂ R 3 × 3 × R 3 is an in variant set of the sy stem dyn amics (20). T o stably embed the system dyn amics ( 2 0) into R 3 × 3 × R 3 , let us consider a function W × R 3 ∋ ( X , Ω) 7→ V ( X , Ω) : = 1 4 k X ⊤ X − I k 2 ∈ R , (21) where W : = { X ∈ R 3 × 3 | det X > 0 } . Then, V − 1 (0) = SO(3) × R 3 and ∇ X V ( X , Ω) = X ( X ⊤ X − I ) , ∇ Ω V ( X , Ω) = 0 . It is easy to sh ow th at V satisfies Assump tion 1, which will actually be proven in th e proo f of The orem 5 . W ith th e help of the func tio n V , th e system dynam ics (20) is exten d ed to the Euclid ean space R 3 × 3 × R 3 as de fin ed in (2) to b e ˙ X = X ˆ Ω − αX ( X T X − I ) , ( 2 2a) ˙ Ω = I − 1 ( I Ω × Ω) + I − 1 u, (22b) where ( X , Ω) ∈ R 3 × 3 × R 3 . T ake a referenc e tr a jectory R ∋ t 7→ ( R 0 ( t ) , Ω 0 ( t )) ∈ SO(3) × R 3 (23) and th e corr espondin g contro l signal R ∈ t 7→ u 0 ( t ) ∈ R 3 such that the trajecto ry obeys th e system dy namics (20). It is trivial to show that R 0 ( t ) satisfies Assumption 2. Define the error trajecto ry as R ∋ t 7→  E ( t ) , Ξ( t )  : =  X ( t ) R − 1 0 ( t ) − I , Ω( t ) − Ω 0 ( t )  ∈ R 3 × 3 × R 3 such that E = 0 , Ξ = 0 en sures that the system dynam ics follows the refe r ence trajectory . The linearized erro r dynamics along the reference state-con tr ol trajectory ( R 0 , Ω 0 , u 0 ) ∈ SO(3) × R 3 × R 3 is ther e fore given by ˙ E = R 0 ˆ Ξ R − 1 0 − 2 α Sy m ( E ) , (24a) ˙ Ξ = I − 1 ( I Ξ × Ω 0 + I Ω 0 × Ξ) + I − 1 δ u, (24b) where Sym ( E ) is the symmetric co m ponen t o f the m atrix E and δ u ( t ) : = u ( t ) − u 0 ( t ) . Now we are in the p osition to split the erro r d ynamics (24 a) into the parallel and the transversal direction. Th e para llel dir ection is given by the Lie a lg ebra so (3) ( the set of 3 × 3 ske w sy mmetric matrices) of the L ie group SO(3) and the transversal d irection is given by the perpen d icular space so (3) ⊥ to th e Lie algebr a so (3) in R 3 × 3 under the Euc lid ean norm, i.e., the set of 3 × 3 symme tric matrices. Consequently , the attitud e er ror d ynamics (24a) is split into the para llel and th e transversal direction , simp lifying (24) to ˙ E ⊥ = − 2 αE ⊥ , (25a) ˙ E k = R 0 ˆ Ξ R − 1 0 , (25b) ˙ Ξ = I − 1 ( I Ξ × Ω 0 + I Ω 0 × Ξ) + I − 1 δ u, (25c) where R ∋ t 7→ E ⊥ ( t ) ∈ so (3) ⊥ and R ∋ t 7→ E k ( t ) ∈ so (3) . Remark 3 . I t is worth noting that the parallel err or dynamics (25b) and the transversal err or dynamics (25a) ar e decoup led. Ther efor e, in the ab sence of a drift vector field, i.e., α = 0 , the initial err or in th e transversal dir ection can not be mitigated and that lead s to a steady-state err or in the transversal dir ection. In other wor ds, for α = 0 , the lin earized err or dynamics (25 a) can not be sta b ilized to zer o. The discre tize d dyn amics of (25), by Euler’ s meth od, is giv en by E ⊥ k +1 = E ⊥ k − 2 h αE ⊥ k , (26a) E k k +1 = E k k + hR 0 ,k ˆ Ξ k R − 1 0 ,k , (26b) Ξ k +1 = Ξ k + h I − 1 ( I Ξ k × Ω 0 ,k + I Ω 0 ,k × Ξ k ) + I − 1 δ u k , (26c) 6 where h is the sam pling time. The following th eorem proves exponential stability of (2 6a): Theorem 5. The transversal err or dyna mics (2 6a) is e xpon e n- tially stab le if 0 < αh < 1 . (27) Pr oof. W e e m ploy Theo rem 3 to establish the stability of the dy namics ( 26a). It is easy to p rove tha t the function V in (21) satisfies V ( X R ) = V ( X ) f or all X ∈ R 3 × 3 , R ∈ SO (3 ) and ∇ 2 V ( I )( X ⊥ , X ⊥ ) = 2  X ⊥ , X ⊥  for all X ⊥ ∈ so (3) ⊥ . Therefo r e the min imum eigenv alue λ min and the maximum eigenv alue λ max of the operato r ∇ 2 V ( I ) ar e equal to 2. Furth er , using the fact that R ⊤ 0 ( t ) R 0 ( t ) = I for all t lead s to β g min = β g max = 1 . Applying Th eorem 3 to the dynamics (2 6a), we obtain (27). Remark 4. In an iden tica l manner , one can pr ove e xponentia l stabilizability of the system dyn a mics (26) by app lying Theo- r em 4. A. Mo del predictive contr ol design In this section, we design a mod el pred ictiv e tracking control of the discrete-time a ttitu de con trol dy namics (2 6). Notice that the transversal dy namics (26a) in ( 26) is deco upled from (26b) and (26c) and expon entially stable (see Theorem 5). Th erefor e , it is advantageous to choose Q E , Q f E in ( 17) which d ecouples the c o st ( 17) along the parallel and the transversal directio n, i.e., k E N k 2 Q f E = k E k N k 2 Q f E + k E ⊥ N k 2 Q f E and k E i k 2 Q E = k E ⊥ i k 2 Q E + k E k i k 2 Q E for all k = 0 , . . . , N − 1 , so that we can ignore the transversal dynam ics as it is not influencing th e optim ization prob lem. Consequently , an N h orizon optimal contro l problem (19) at a g iv en time instant k for the system dyn amics (26) with the perfor mance objective (17) and con stra in ts (18) is given by minimize { δu k + i | k } N − 1 i =0 J  E k k : k + N | k , Ξ k : k + N | k , δ u k : k + N | k  subject to                                       E k k + i +1 | k = E k k + i | k + hR 0 ,k + i | k ˆ Ξ k + i | k R − 1 0 ,k + i | k Ξ k + i +1 | k = Ξ k + i | k + h I − 1 ( I Ξ k + i | k × Ω 0 ,k + i | k ) + h I − 1 ( I Ω 0 ,k + i | k × Ξ k + i | k ) + I − 1 δ u k + i | k δ u k + i | k ∈ U k + i for all i = 0 , . . . , N − 1 , ( E k k + i | k ∈ X E k + i Ξ k + i | k ∈ X Ξ k + i for all i = 1 , . . . , N ,  E k k | k , Ξ k | k  =  E k k , Ξ k  (28) where  E k k , Ξ k  is fixed, and X E k + i and X Ξ k + i are admissible sets for E k k + i | k and Ξ k + i | k , respectively . Th e quadra tic pro - gram (28) can be solved in MA TLAB using an optimization modeling too lb ox Y ALMIP [1 8]. A detailed exposition of computatio nal com plexity of real- tim e MPC exists in [1 9]– [21]. T h e op timal co n trol pro blem (28) is so lved at each time instant k and the contro l u : = u 0 + δ u, where δ u : = δ u k | k for [ k h, ( k + 1) h [ , is applied to th e system as shown in Figure 1 . ZOH Attitude dy namics Disturbances Sampler MPC ( E k , Ξ k ) δ u + u ( E , Ξ) ( X 0 , Ω 0 ) − ( X, Ω) u 0 + + δ u k Fig. 1: T he sampled -data closed-loo p system with MPC. I V . S I M U L A T I O N R E S U LT S W e simulate o ur MPC in a sampled- d ata system as shown in Figure 1. T h e finite d imensional quadr atic pro grammin g problem (28) is solved at each time instan t k to calculate a feedback control law . The moment of in ertia matr ix of th e rigid b ody in (28) is I = diag (4 . 250 , 4 . 3 37 , 3 . 66 4) , which was taken fr o m a satellite f r om the Euro p ean Student Earth Orbiter (ESEO) [22]. L et us track a refere n ce tr ajectory R ∋ t 7→ ( R 0 ( t ) , Ω 0 ( t )) ∈ SO(3) × R 3 , where R 0 ( t ) : = exp ( t ˆ e 1 ) exp ( t ˆ e 2 ) exp ( t ˆ e 3 ) with e i as the un it vector along the i th axis, and Ω 0 ( t ) : =  1 + cos t, sin t − sin t cos t, cos t + sin 2 t  ⊤ , using the MPC c o ntrol law with the co rrespond ing refer ence control signa l R ∋ t 7→ u 0 ( t ) = I ˙ Ω 0 ( t ) − ( I Ω 0 ( t )) × Ω 0 ( t ) ∈ R 3 . The initial d ata and par ameters considered for the optimiza- tion (2 8) are the f o llowing: • E k 0 =  R 0 (0 . 2) − R 0 (0)  k , Ξ 0 = (0 , 0 , 0) ⊤ , • Q E = 100 I , Q Ξ = 10 I , Q δu = 0 . 01 I , • Q f E = 100 I , Q f Ξ = 10 I , • sampling time: h = 0 . 2 sec , • MPC time horizon : N = 4 . The MPC controller takes the tr a cking erro r measuremen ts  E k , Ξ k  for comp uting the fe e d back co ntrol δ u k at the each iteration k . The se track ing error measure m ents are calcu lated as a d ifference o f the re f erence states  R 0 ,k , Ω 0 ,k  and th e states ob tained fro m the ODE simulation (we have used the MA TLAB integrator, ode45, with the option s, Re l T ol = AbsT ol = 10 − 6 ) of th e rigid body dynam ics (20) . In turn, the ODE simulation is d r iv en by the zer o-ord er ho ld con trol actions gener ated by the MPC contr oller , and that forms the closed-loo p M PC sy stem; see Figu re 1. W e simulate th ree different scenarios as follows: 7 A. Ca se 1: Loose constraint and no n oise In the first case study , we con sid e r the f ollowing state a n d control con straints (18): X E k = so (3) , X Ξ k = R 3 , U k = U − u 0 ,k for ea c h k , where U : = { y ∈ R 3 | − 10 ≤ y i ≤ 10 for i = 1 , 2 , 3 } . The closed-loo p system with the d esigned MPC shows a successful tracking perfo r mance as the erro r trajectory ( E ( t ) , Ξ( t )) tend s to zero quic k ly (see Figure 2 ), an d the optima l co ntrol pr ofile obeys th e con trol constraints as shown in Figure 3. It is worth noting tha t the angular velocity erro r Ξ shown in Figu re 2 increases initially f r om z e r o in or der to m itigate the initial orientation er ror E . 0 2 4 6 8 10 0 0.2 0.4 0.6 0 2 4 6 8 10 0 0.2 0.4 0.6 0.8 P S f r a g r e p l a c e m e n t s T ime [ sec ] T ime [ sec ] Fig. 2: The tracking er rors fo r Case 1. 0 2 4 6 8 10 -10 0 10 0 2 4 6 8 10 -10 0 10 0 2 4 6 8 10 -10 0 10 P S f r a g r e p l a c e m e n t s T ime [ sec ] T ime [ sec ] T ime [ sec ] T i m e [ s e c ] Fig. 3: Zero-o rder hold contro l for the sampled data system with M PC for Case 1 . B. Ca se 2: T i ght contr ol con straint The second case we study is the one with th e following tight con straints: X E k = so (3) , X Ξ k = R 3 , U k = U − u 0 ,k for ea c h k , where U : = { y ∈ R 3 | − 6 ≤ y i ≤ 6 for i = 1 , 2 , 3 } . T he closed-loo p system with th e designed MPC considering a tigh t control bound shows a compro mised tracking performan ce. As the co n trol trajectory h its the c o ntrol bo unds (see Figu re 5), the error trajectory ( E ( t ) , Ξ( t )) deviates fro m zer o , as sh own in Figur e 4. 0 2 4 6 8 10 0 0.2 0.4 0.6 0 2 4 6 8 10 0 0.2 0.4 0.6 0.8 P S f r a g r e p l a c e m e n t s T ime [ sec ] T ime [ sec ] Fig. 4: The tracking er rors fo r Case 2. 0 2 4 6 8 10 -10 0 10 0 2 4 6 8 10 -10 0 10 0 2 4 6 8 10 -10 0 10 P S f r a g r e p l a c e m e n t s T ime [ sec ] T ime [ sec ] T ime [ sec ] Fig. 5: Zero-o rder hold contro l for the sampled data system with MPC for Case 2 . C. Case 3 : Noisy measur ements In this case stud y , the state and control constrain ts are con - sidered as in Case 1 in Section IV -A; howev er , a m e a surement noise in the tracking er r or  E k k + i +1 | k , Ξ k + i +1 | k  ∈ so (3) × R 3 for i = 1 , . . . , N , is realized by the indepen dent and identically distributed (i.i. d.) random variables, i.e., w ∼ N (0 , σ 2 w ) , wher e σ w = 0 . 03 . Due to the noisy state m e asurements, the erro r trajecto ry ( E ( t ) , Ξ( t )) fluctua te s aro und zero instead of stabilizin g at zero; see Figur e 6. However , th e tr a cking per forman ce of the 8 closed lo op system is similar to Case 1, as shown in Figure 6 and Figu re 7. 0 2 4 6 8 10 0 0.2 0.4 0.6 0 2 4 6 8 10 0 0.2 0.4 0.6 0.8 P S f r a g r e p l a c e m e n t s T ime [ sec ] T ime [ sec ] Fig. 6: The tracking er rors fo r Case 3. 0 2 4 6 8 10 -10 0 10 0 2 4 6 8 10 -10 0 10 0 2 4 6 8 10 -10 0 10 P S f r a g r e p l a c e m e n t s T ime [ sec ] T ime [ sec ] T ime [ sec ] Fig. 7: Zero-o rder hold contro l for the sampled data system with M PC for Case 3 . V . C O N C L U S I O N In this pap er , we hav e presented a tech n ique to design model predictive tra c king co ntrollers for con trol systems ev olving on manifold s in Euclid ean spaces. W e h av e app lied the propo sed technique to the systems on matrix Lie groups and demon- strated its potency by designing a linear MPC law for th e rigid body attitu de d y namics. Our ap proach simplifies MPC design for contr ol systems on man ifolds. This development could be quite useful f o r control enginee r s in dealing with nonlinear mechanical co ntrol system applicatio n s in p ractice. R E F E R E N C E S [1] F . Borrelli, A. Bemporad, and M. Morari, Predict ive contr ol for linear and hybrid systems . Cambridge Univ ersity Press, 2017. [2] ——, Predic tive contr ol for linea r and hybrid systems . Cambridge Uni versi ty Press, 2017. [3] J . B. Rawlings, D. Q. Mayne, and M. Diehl, Model Pre dictiv e Contr ol: Theory , Computation, and Design . Nob Hill Publishing, 2017. [4] S . H. Lee and C. C. Chung, “Multile v el approximate model predicti v e control and its applicati on to autonomous vehic le acti ve steering, ” in IEEE 52nd Annual Confer ence on Decision and Contr ol (CDC) . IE EE, 2013, pp. 5746–5751. [5] H . Li, W . Y an, and Y . Shi, “Continu ous-time model predicti ve control of under-actu ated spacecraft with bounded control torques, ” Automatica , vol. 75, pp. 144–153, 2017. [6] B. Stellato, T . Geyer , and P . J. Goulart , “High-spee d finite control set model predict i ve control for power electroni cs, ” IEEE T r ansactions on P o wer Electr oni cs , vol. 32, no. 5, pp. 4007–4020, May 2017. [7] D . E. Chang, “On contro ller design for systems on manifolds in Euclidea n space, ” International J ournal of Robust and Nonlin ear Contr ol , vol. 28, no. 16, pp. 4981–4998, 2018. [Online]. A v ailabl e: https:/ /doi.org/ 10.1002/rnc.429 4 [8] D . E. Chang, F . Jim ´ enez, and M. Perlmutter , “Feedback integrat ors, ” J . Nonlinear Sci. , vol. 26, no. 6, pp. 1693–1721, 2016. [Online]. A vaila ble: https://do i.org/10 .1007/s00332- 016- 9316- 7 [9] U . V . Kalabi ´ c, R. Gupta, S. Di Cairano, A. M. Bloch, and I. V . Ko lmanovsk y , “MPC on manifol ds with an appl icati on to the cont rol of spacecr aft attitude on SO(3) , ” Automatica , vo l. 76, pp. 293–300, 2017. [Online]. A v ai lable: https:/ /doi.org/ 10.1016/j.automat ica.2016.10.022 [10] A. Kny aze v and A. Malyshev , “Continu ation model predicti v e control on smooth m anifolds, ” IF AC -P ape rsOnLine , vol . 48, no. 25, pp. 126 – 131, 2015, 16th IF A C W orkshop on Control Applica tions of Optimization CA O2015. [Online]. A v ail able: http:/ /www .scie ncedirec t.com/science/article/pii/S2405896315023290 [11] D. Nesic, A. T e el, and P . Kokot ovic, “On the design of a controll er based on the discr ete-time approximation of the nonlinea r plant model, ” in Proce edings of the 1999: American Contr ol Conferen ce , vol. 5. IE EE, 1999, pp. 3474–3478. [12] D. Nesic and A. R. T ee l, “ A framew ork for s tabili zation of nonlinea r sampled-da ta systems based on their approximate discrete -time m odels, ” IEEE T ransact ions on Automat ic Contr ol , vol. 49, no. 7, pp. 1103–1122, 2004. [13] R. Abraham, J. E. Marsden, and T . Ratiu, Manifol ds, tensor analysis, and applicati ons . S pringer Science & Business Media, 2012, vol. 75. [14] S. Sastry , Nonlinear Systems: Analysis, Stabil ity , and Contr ol . Springer Science & Business Media, 2013, vol. 10. [15] A. M. Bloch, Nonholonomi c mechanic s and contr ol , 2nd ed., ser . Interdisc iplinar y Applie d Mathemati cs. Springe r , New Y ork, 2015, vol. 24, with the collaborat ion of J. Bailieul, P . E. Crouch, J. E. Marsden and D. Z enko v , Wi th scientific input from P . S. Krishnaprasad and R. M. Murray . [Online]. A vaila ble: https://doi.or g/10.1007/ 978- 1- 4939- 3017- 3 [16] J. B. Rawli ngs and K. R. Muske, “The stabilit y of constra ined recedi ng horizon control, ” IEEE T ransac tions on Automatic Contr ol , vol. 38, no. 10, pp. 1512–151 6, 1993. [Online]. A vaila ble: https:/ /doi.org/ 10.1109/9.241565 [17] P . Falcone , F . Borrell i, H. E . T seng, J. Asgari, and D. Hrov at, “Linear time-v aryin g model predicti v e control and its app licati on to activ e steering systems: Stability analysis and expe rimental v alidat ion, ” Inter- national Journal of Robust and Nonlinear Contr ol , vol. 18, no. 8, pp. 862–875, 2008. [18] J. L ¨ ofbe rg, “Y ALMIP: A toolbox for m odeling and optimiz ation in MA TLAB, ” in Pr oceedi ngs of the CACSD Confer ence , vol. 3. T aipei, T a iwan, 2004. [19] S. Richte r , C. N. Jones, and M. Morari, “Computati onal complexit y certific ation for real-ti me MPC with input constrai nts based on the fast gradien t m ethod, ” IEEE T ransact ions on Automat ic Contr ol , vol. 57, no. 6, pp. 1391–1403, 2012. [20] M. Rubagot ti, P . Patrinos, and A. Bemporad, “Stabiliz ing embedded MPC with computati onal comple xity guarant ees, ” in 2013 Eur opean Contr ol Confer ence (ECC) . IEEE, 2013, pp. 3065–3070 . [21] V . Nedelcu, I. Necoara , and Q . Tran-Dinh, “Computati onal complexity of ine xact gradie nt augmented Lagrangi an methods: appli cation to constrai ned MPC, ” SIAM Journal on Contr ol and Optimization , vol. 52, no. 5, pp. 3109–3134 , 2014. [Online]. A v ail able: https:/ /doi.org/ 10.1137/12089 7547 [22] Ø. Hegrenæs, J. T . Gravdahl , and P . Tøndel, “Space craft attitude control using expl icit model predic ti ve control, ” Automatic a , vol. 41, no. 12, pp. 2107–211 4, 2005. [Online]. A vaila ble: https:/ /doi.org/ 10.1016/j.aut omatica.2005.06.015